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In  het  kader  van  het  project  ‘Feasability  study  of  an  IM-Ignitor’  is  onderzoek 
gedaan  naar  de  ontwikkeling  van  een  nieuwe  ontsteker  voor  een  raketsysteem.  Het 
doel  van  dit  onderzoek  was  een  ontsteker  te  ontwikkelen  met  een  lagere  ontsteek- 
temperatuur  zodat  dat  de  kans  voor  raketsystemen  op  een  IM-typeclassificatie 
vergroot  zou  worden. 

In  het  bovengenoemde  project  zijn  veel  uiteenlopende  aspecten  onderzocht.  Een 
groot  deel  daarvan  betrof  het  onderwerp  ‘Cook-off’.  Door  enerzijds  de  tijdslimiet 
en  anderzijds  budgettaire  redenen  konden  niet  alle  aspecten  op  een  gewenste  ma- 
nier  afgehandeld  worden.  Omdat  een  aantal  zaken  interessant  is  voor  het  project 
‘Thermische  Initiatie’  (TI  onder  opdrachtnummer  A95KL408)  is  in  het  kader 
hiervan  in  het  vorige  en  het  huidige  jaar  ruimte  gecreeerd  om  een  aantal  van  die 
aspecten  te  evalueren  en  af  te  ronden.  De  samenvatting  van  de  resultaten  en  de 
evaluatie  van  het  bovengenoemde  ‘Feasibility’ -project  zijn  reeds  gerapporteerd  [3]. 
De  interessante  ontwikkelingen  op  het  gebied  van  computermodeliering  en  de 
resultaten  van  een  aantal  temperatuur-,  reactiekinetiek-  en  drukberekeningen  en 
nieuwe  berekeningen  met  de  gegevens  uit  het  feasibility-project  worden  in  dit 
rapport  beschreven. 

Er  is  een  aantal  computermodellen  ontwikkeld  waarmee,  ofwel  eenvoudige  ther- 
misch-chemische  berekeningen  ofwel  meer  complexe  gecombineerde  thermisch- 
chemische  drukberekeningen  uitgevoerd  kunnen  worden.  Naast  gegevens  uit  de 
literatuur  zijn  thermische-  en  reactiesnelheidsparameters  bepaald  die  nodig  zijn 
voor  deze  modellen.  Met  de  ontwikkelde  computerprogramma’s  zijn  vervolgens 
TNO  Prins  Maurits  Laboratorium  (TNO-PML)  Cook-off  buis-,  raketsysteem-  en 
slow  Cook-off-berekeningen  uitgevoerd  en  vergeleken  met  de  experimentele 
waarden. 

Met  de  resultaten  van  de  temperatuurprofielberekeningen  van  het  raketsysteem  is 
een  beter  inzicht  verkregen  in  het  warmtetransport  door  het  systeem.  Dit  kan  bij  de 
ontwikkeling  van  nieuwe  systemen  een  bijdrage  leveren  aan  het  ontwerp.  Een 
betere  inschatting  van  een  mogelijke  responsie  van  een  ratketsysteem  op  een  ther¬ 
mische  stimulus  kan  gemaakt  worden  indien  naast  deze  berekeningen  de  resultaten 
van  de  verschillende  Cook-off  experimenten  toegevoegd  worden. 
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Het  blijkt  dat  de  reactiesnelheidsvergelijkingen  die  bepaald  zijn  voor  de  stuwstof- 
fen  HTPB/AP  en  PPG/AP/AN  niet  overeenkomen  met  de  werkelijkheid.  Hiervoor 
wordt  een  aantal  redenen  genoemd  waaruit  blijkt  dat  de  bepaling  van  de  reactie- 
snelheidvergelijking  van  energetische  materiaal  complexer  is  dan  in  eerste  instan- 
tie  aangenomen  werd. 

Met  het  gebruik  van  de  computermodellen  met  gegevens  uit  de  literatuur,  de 
Tarver-McGuire  reactiekinetiek  en  gecombineerde  drukberekeningen  zijn  goede 
resultaten  bereikt.  Ook  uit  eerdere  berekeningen  is  gebleken  dat  aan  de  hand  van 
dit  model  met  nulde  orde  kinetiek  al  goede  resultaten  behaald  konden  worden  [26]. 
De  temperatuurprofielberekeningen  geven  een  goede  indruk  van  de  temperatuur- 
gradienten  van  het  testitem  die  zeer  goed  overeenkomen  met  de  experimentele 
waarden.  Met  behulp  van  de  drukberekeningen  wordt  een  stopcriterium  voor  de 
simulaties  verkregen  die  beter  met  de  werkelijkheid  overeenkomen  dan  slechts 
zuiver  thermische-chemische  berekeningen. 
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Samenvatting 


In  het  kader  van  het  project  ‘Feasibility  study  of  an  IM-Ignitor’  is  onderzoek 
gedaan  naar  de  ontwikkeling  van  een  nieuwe  ontsteker  voor  een  raketsysteem.  In 
dit  project  zijn  veel  uiteenlopende  aspecten  onderzocht.  Een  groot  deel  daarvan 
betrof  het  onderwerp  Cook-off  en  is  daarom  interessant  voor  het  project 
‘Thermische  Initiatie’.  De  interessante  ontwikkelingen  op  het  gebied  van  compu- 
termodellering  en  de  resultaten  van  een  aantal  temperatuur-,  reactiekinetiek-  en 
drukberekeningen  en  nieuwe  berekeningen  met  de  gegevens  uit  het  feasibility 
project  worden  in  dit  rapport  uitgebreid  beschreven.  Er  is  een  aantal  computermo- 
dellen  ontwikkeld  waarmee,  ofwel  eenvoudige,  thermisch-chemische  berekeningen 
ofwel  meer  complexe  gecombineerde  thermisch-chemische-drukberekeningen 
uitgevoerd  kunnen  worden.  Met  de  ontwikkelde  computerprogramma’s  zijn  ver- 
volgens  TNO  Prins  Maurits  Laboratorium  (TNO-PML)  Cook-off  buis-,  raketsys¬ 
teem-  en  slow  Cook-off-berekeningen  uitgevoerd  en  vergeleken  met  de  experimen- 
tele  waarden.  Met  de  resultaten  van  de  temperatuurprofielberekeningen  van  het 
raketsysteem  is  een  beter  inzicht  verkregen  in  het  warmtetransport  door  het  sys- 
teem.  Dit  kan  bij  de  ontwikkeling  van  nieuwe  systemen  een  bijdrage  leveren  aan 
het  ontwerp.  Verder  is  gebleken  dat  de  bepaling  van  de  reactiesnelheidvergelijking 
van  energetische  materiaal  complexer  is  dan  in  eerste  instantie  aangenomen  werd. 
Met  het  gebruik  van  de  computermodellen  met  gegevens  uit  de  literatuur,  de 
Tarver-McGuire  reactiekinetiek  en  gecombineerde  drukberekeningen  zijn  goede 
resultaten  bereikt.  De  temperatuurprofielberekeningen  geven  een  goede  indruk  van 
de  temperatuurgradienten  van  het  testitem  die  zeer  goed  overeenkomen  met  de 
experimentele  waarden.  Met  behulp  van  de  drukberekeningen  wordt  een  stopcrite- 
rium  voor  de  simulaties  verkregen  die  beter  met  de  werkelijkheid  overeenkomen 
dan  slechts  zuiver  thermische-chemische  berekeningen. 
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1  Inleiding 


In  het  kader  van  het  project  ‘Feasability  study  of  an  IM-Ignitor’  is  onderzoek 
gedaan  naar  de  ontwikkeling  van  een  nieuwe  ontsteker  voor  een  raketsysteem.  Het 
doel  van  dit  onderzoek  was  een  ontsteker  te  ontwikkelen  met  een  lagere  ontsteek- 
temperatuur  zodat  dat  de  kans  voor  raketsystemen  op  een  IM-typeclassificatie 
vergroot  zou  worden. 

In  het  bovengenoemde  project  zijn  veel  uiteenlopende  aspecten  onderzocht.  Een 
groot  deel  daarvan  betrof  het  onderwerp  ‘Cook-off.  Door  enerzijds  de  tijdslimiet 
en  anderzijds  budgetaire  redenen  konden  niet  alle  aspecten  op  een  gewenste  ma- 
nier  afgehandeld  worden.  Omdat  een  aantal  zaken  interessant  is  voor  het  project 
‘Thermische  Initiatie’(TI  onder  opdrachtnr  A95KL408)  is  in  het  kader  van  in  het 
vorige  en  het  huidige  jaar  ruimte  gecreeerd  om  een  aantal  van  die  aspecten  evalue- 
ren  en  af  te  ronden.  De  samenvatting  van  de  resultaten  en  de  evaluatie  van  het 
bovengenoemde  ‘Feasibility’-project  zijn  reeds  gerapporteerd  [3].  De  interessante 
ontwikkelingen  op  het  gebied  van  computermodellering  en  de  resultaten  van  een 
aantal  temperatuur-,  reactiekinetiek-  en  drukberekeningen  en  nieuwe  berekeningen 
met  de  gegevens  uit  het  ‘Feasibility’-project  worden  in  dit  rapport  beschreven. 

Algemeen  is  bekend,  dat  onder  meer  de  slow  Cook-off  test  problemen  oplevert 
voor  een  IM-classificering  van  munitie-items.  Bij  langzame  opwarming  bevindt  de 
gehele  stof  zich  op  een  verhoogde  temperatuur,  zodat  de  mechanische  eigenschap- 
pen  snel  verslechteren  en  de  poreusheid  toeneemt.  In  het  geval  van  een  slow  Cook¬ 
off  reactie  zal  de  thermische  wegloopreactie  veelal  vanuit  het  midden  van  de 
explosive  stof  starten  waardoor  een  extra  opsluiting  van  de  stof  optreedt.  Door  de 
combinatie  van  de  verslechterde  mechanische  eigenschappen,  verhoogde  poreusi- 
teit  en  de  interne  opsluiting  van  de  explosieve  stof  in  combinatie  met  een  opslui¬ 
ting  van  het  omhullende  kan  een  verbrandingsreactie  overgaan  in  een  detonatie 
(DDT-overgang). 

In  het  ‘Feasibility’-project  is  gezocht  naar  de  mogelijkheid  een  raketsysteem 
minder-gevoelig  te  maken  voor  een  (slow)  Cook-off  test.  Om  dit  doel  te  bereiken 
werd,  ten  tijde  van  het  project,  uitgegaan  van  de  filosofie  dat,  indien  een  ontsteker 
bij  een  lagere  temperatuur  dan  de  hoofdstuwstof  ontsteekt  (ontleedt),  de  raket  op 
een  normale  wijze  kan  functioneren  mits  de  hoofdlading  bij  deze  ontsteektempera- 
tuur  nog  een  normaal  brandgedrag  vertoont.  Verder  is  bij  deze  filosofie  de  aanna- 
me  gedaan  dat  de  hoofdlading  bestaat  uit  een  HTPB/AP-stuwstof  omgeven  door 
een  metalen  casing.  De  keuze  voor  de  ontstekerstuwstof  is  gevallen  op  een 
PPG/AP/AN  stuwstof  die,  ten  gevolge  van  de  toevoeging  van  AN,  bij  een  lagere 
temperatuur  zou  moeten  ontleden  en  ontsteken. 

Om  in  een  Cook-off  test  de  ontsteker  vroeger  te  laten  ontsteken  dan  de  hoofdlading 
moet  de  warmtestroom  in  de  raket  zodanig  zijn  dat  de  ontsteker  eerder  op  ontsteek- 
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temperatuur  is  dan  de  hoofdlading.  Daar  de  ontsteker  zich  in  een  raket  altijd  in  het 
centrum  bevindt,  wordt  niet  vanzelf  aan  deze  eis  voldaan.  Om  beter  inzicht  te 
krijgen  in  de  temperatuurprofielen  in  een  raket,  zijn  computersimulaties  uitgevoerd 
om  deze  warmteprofielen  bij  langzame  en  intermediate  opwarmsnelheid  te  bereke- 
nen. 

Er  is  tijdens  het  ‘Feasibility  project’  voor  een  stap  voor  stap  benadering  van  de 
oplossing  voor  dit  complexe  probleem  gekozen.  Hierdoor  zou  bij  uitloop  of  niet 
bruikbare  modelparameters  van  andere  werkpakketten  met  de  reeds  gegenereerde 
resultaten  toch  een  zo  goed  mogelijke  conclusie  getrokken  kunnen  worden.  In  het 
vervolg  van  het  rapport  zal  over  ‘een  eerste  fase’  gesproken  worden  indien  het 
gegevens  uit  de  ‘Feasibility’-studie  aangaat,  en  over  ‘een  tweede  fase’  indien  het 
vervolgwerkzaamheden  betreft  die  in  het  kader  van  het  project  ‘Thermische  Initia¬ 
te  ’  uitgevoerd  zijn. 

Als  eerste  stap,  in  fase  een,  werd  een  model  van  de  TNO  Prins  Maurits  Laborato- 
rium  (TNO-PML)  Cook-off  buis  met  een  HMX-PBX  [4]  doorgerekend  en  vergele- 
ken  met  experimentele  resultaten.  In  een  tweede  stap  werd  na  een  uitgebreide 
analyse  van  verschillende  raketsystemen  een  rekengrid  van  het  raketsysteem 
gemaakt.  Met  dit  model  werden  in  een  derde  stap  thermische  profielen  berekend. 
Tevens  werden  de  opwarmsnelheid  en  de  geleidingscoefficient  van  het  materiaal, 
dat  de  warmtestroom  in  de  richting  van  de  ontsteker  bepaald,  gevarieerd. 

Eveneens  is  met  behulp  van  thermogravimetrie  de  reactiekinetiek  van  de  kruiten 
bepaald  en  zijn  TNO-PML  Cook-off  experimenten  met  de  energetische  materialen 
uitgevoerd.  Ook  is  een  aantal  thermische  parameters  zoals  de  soortelijke  warmte 
en  de  warmtegeleidingscoefficient  bepaald. 

Met  behulp  van  de  ontwikkelde  programmatuur  en  gegevens  zijn,  in  de  tweede 
fase,  in  het  kader  van  het  project  ‘Thermische  Initiatie’  berekeningen  uitgevoerd. 
Hieruit  zou  moeten  blijken  of  de  berekeningen  met  het  model  slechts  een  globale 
indruk  van  de  temperatuurprofielen  zouden  kunnen  verschaffen  dan  wel  een  goede 
schatting  van  de  druk,  explosietijd  en  -temperatuur.  Met  de  reactiekinetiek  zullen 
tenslotte  twee  berekeningen  uitgevoerd  worden  aan  het  raketmodel. 

In  hoofdstuk  2  wordt  de  theorie  uiteengezet  omtrent  het  drukberekeningsmodel,  in 
hoofdstuk  3  wordt  de  implementatie  van  het  model  in  de  computercode  ABAQUS 
[2]  beschreven,  worden  de  gridmodellen  toegelicht  en  de  verrichtte  werkzaamhe- 
den  besproken.  In  hoofdstuk  4  worden  de  resultaten  van  fase  1,  het  ‘Feasibilty’- 
project,  gegeven,  gevolgd  door  de  resultaten  van  fase  2  in  hoofdstuk  5. 

Hoofdstuk  6  bevat  de  discussie  gevolgd  door  de  conclusies  in  hoofdstuk  7. 
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2  Theorie 


2.1  Het  drukmodel 

Het  gebruikte  model  is  een  gedeelte  van  een  model  zoals  beschreven  staat  in  onder 
andere  ‘A  Constitutive  mechanical  model  for  energetic  materials’  van 
M.L.  Hobbss,  M.R.  Baer  en  R.J.  Gross  [6].  Het  model  gaat  uit  van  een  ‘bubble’ 
(bellen)  mechanisme,  hetgeen  inhoudt  dat  bij  de  decompositie  van  veel  energeti- 
sche  materialen  holten  in  de  explosieve  stof  ontstaan  ten  gevolge  van  de  omzetting 
van  vaste  stof  in  gas.  Aan  de  hand  van  een  beginporeusiteit,  de  hoeveelheid  omge- 
zette  stof  en  de  vrijgekomen  hoeveelheid  gas,  kan  een  schatting  gemaakt  worden 
van  de  lokale  druk  in  deze  holten.  Daar  er  van  een  poreusiteit  van  de  explosieve 
stof  uitgegaan  kan  worden,  vindt  diffusie  van  het  gas  plaats  waardoor  er  sprake  is 
van  een  globale  drukstijging  in  de  buis. 

Ten  gevolge  van  een  beginporeusiteit  bevat  een  energetisch  materiaal  een  aantal 
microbellen  (bubbles)  per  eenheid  van  ruimte  (cel  of  computerelement).  Door  de 
ontledingsreactie  van  de  explosieve  stof  zullen  deze  microbellen  groeien  in  deze 
poreusiteiten,  ook  wel  defecten  genoemd.  In  het  model  hebben  de  bellen  een 
exteme  diameter  b ;  een  interne  diameter  a  en  een  wanddikte  van  (b-a).  In  dit 
eenvoudige  model  stelt  b  de  gemiddelde  afstand  van  deze  defecten  voor  en  is  a  de 
grootte  van  het  defect.  De  inwendige  bol  bevat  gas  met  een  initiele  dichtheid  pg° 
en  wordt  gevuld  met  decompositiegassen.  De  begindichtheid  van  de  explosieve 
stof  is  pc°.  The  initiele  binnenstraal  a0  en  buitenstraal  ba  van  de  bellen  kunnen 
bepaald  worden  aan  de  hand  van  initiele  gasvolumefractie  <|>0  en  de  bellendichtheid 
N0  met  de  volgende  formules 

«o=fco<l)o1/3  (1) 

en 

b0=(3/ 4kN0)113  (2) 

met  4>o  de  initiele  gasvolumefractie  en  No  de  bellendichtheid  in  het  materiaal. 

Omdat  slechts  een  gedeelte  van  het  model  gebruikt  wordt  om  de  inwendige  druk 
ten  gevolge  van  gasontwikkeling  te  berekenen,  worden  van  de  zes  algebraische 
hoofdvergelijkingen  er  drie  uit  het  beschreven  model  [6]  gebruikt.  Een  vergelijking 
voor  het  behoud  van  de  gasmassa  (3),  behoud  van  totale  massa  (4)  en  de  gasfase- 
toestandsvergelijkingen  (5). 
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De  vergelijking  voor  behoud  van  de  gasmassa  luidt  als  volgt: 


P* 


P°<t>o¥ 

<t> 


(3) 


Hierin  zijn  pg  en  pg°  respectievelijk  de  dichtheid  en  initiele  dichtheid  van  het  gas, 
y  de  genormaliseerde  gasmassa.  Met  de  aanname  b=bo,  omdat  aangenomen  wordt 
dat  de  buis  nauwelijks  uitzet  bij  druktoename,  volgt: 

<}>  =(1  -tyo)Fe  +  <))o  (3a) 

met  Fe  de  effectieve  waarde  van  de  verhoudingen  van  de  massa  van  de  gasproduc- 
tie  en  de  initiele  vaste  stofmassa  in  de  eenheidscel.  Deze  waarde  is  eenvoudig  te 
bepalen  door  de  berekende  massafractie  van  de  decompositiegassen  F  te  verme- 
ningvuldigen  met  een  dimensieloze  factor  (1-Xc).  De  massafractie  F  wordt  door  de 
reactiesnelheidsvergelijkingen  bepaald.  Op  de  parameter  Xc  wordt  in  de  gasfase- 
vergelijking  nog  dieper  ingegaan.  De  eerdergenoemde  parameter  <t>o  kan  eenvoudig 
bepaald  worden  met  de  vergelijking: 

<|)0  =1 — ^  (3b) 

Pth 


Fe  ziet  er  dan  als  volgt  uit. 

Fe  =  F(l-Xc)  (3c) 

De  materiaalparameter  A.  wordt  als  volgt  gedefmieerd. 

D°  1 

X=%(-!--l)  +  l  (3d) 

P°  0O 


De  vergelijking  voor  het  totale  massabehoud  ziet  er  als  volgt  uit. 

pc=ps(--l)/(i-l)  W 

*  \jr  4> 

met  pc  de  dichtheid  van  de  vaste  stof.  Uiteindelijk  kan  met  behulp  van  de  volgende 
vergelijking,  de  gasfasetoestandsvergelijking,  de  druk  Pg  in  een  eenheidscel  bere- 
kend  worden. 

Pg  =  zPgRT  /  Mw  (5) 

met  R  de  universele  gasconstante,  T  de  temperatuur  in  de  eenheidscel  en  Mw  is  het 
gemiddelde  moleculaire  gewicht  in  de  gasfase.  De  BKW-compressibiliteitsfactor  z 
is  gedefineerd  als  volgt. 


Z  =  1  +  xe 


(5a) 
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en  de  BKW-parameter  x  is  gedefinieerd  als  volgt. 

x  =  pgK^yiki/(Mw(T  +  0f)  (5b) 

i 

Voor  de  gasfasetoestandsvergelijking  (4)  is  in  het  model  gebruikgemaakt  van  de 
zogenaamde  BKW-toestandsvergelijking.  Deze  voldoet  goed  bij  hoge  gasdrukken 
en  is  gedefineerd  voor  hoge  gasdichtheden.  Met  behulp  van  de  TIGER-code  [7] 
(Cowperthwaite  and  Zwisler,  1973)  en  deze  BKW-toestandsvergelijking  (Hobbs 
and  Baer,  1993)  [8]  kan  het  gemiddelde  gasfase  moleculaire  gewicht  Mw,  het 
covolume  £y;kj  en  massafractie  van  de  eindtoestand  van  de  vloeistof-  en  vaste 
stoffase  van  ontstekingsproducten  Xc,  bepaald  worden.  Het  blijkt  dat,  onder  andere 
voor  de  energetische  materialen  NC,  HMX,  RDX,  TATB  en  TNT  voor  hogere 
drukken  en  temperaturen  waarbij  de  compressibiliteit  z  groter  is  dan  1,  die  eind- 
producten  relatief  constant  zijn  en  daarmee  Xc.  Dit  houdt  tevens  in  dat  het  gemid¬ 
delde  gasfasemoleculaire  gewicht  Mw  en  het  covolume  Eyjkj  vrij  constant  zijn  bij 
hoge  drukken  en  temperaturen.  Met  andere  woorden,  een  waarde  kan  voor  de 
totale  berekening  gebruikt  worden  onder  de  condities  die  verwacht  worden  bij  een 
voorinitiatiefase  in  een  Cook-off  proces. 

2.1.2  Schatting  van  de  globale  druk  in  de  buis 

Nadat  via  bovenstaande  vergelijkingen  de  druk  per  eenheidscel  bepaald  is,  kan  een 
schatting  gemaakt  worden  voor  de  globale  druk  in  de  buis.  Hiervoor  is  slechts 
gebruikgemaakt  van  een  eenvoudige  theorie  om  een  schatting  te  krijgen  van  de 
gemiddelde  druk.  Uitgegaan  wordt  van  n  eenheidscellen  met  elk  een  eigen  tempe- 
ratuur  T„,  een  volume  Vn,  een  druk  Pn  en  een  gasvolumeffactie  <(>n.  In  een  tijdsin- 
terval  zal  ten  gevolge  van  de  decompositie  de  druk  in  een  eenheidscel  veranderen 
omdat  er  vaste  stof  in  gas  omgezet  wordt.  Hierdoor  zal  ook  de  gasvolumefractie 
per  cel  veranderen.  Dit  zal  ten  gevolge  van  de  temperatuurverschillen  in  de  ene  cel 
anders  zijn  dan  de  andere.  Met  behulp  van  de  volgende  formules  wordt  een  schat¬ 
ting  gemaakt  van  de  einddruk  Pe.  In  een  eenheidscel  nemen  we  aan  dat  de  ideale 
gaswet  geldt: 

PnV„=n„RTn  (6) 

Nadat  diffusie  van  het  gas  plaatsgevonden  heeft,  moet  gelden  dat  het  aantal  molen 
gas  voor  en  na  gelijk  is: 

ni  +  n2  +  n3  + +nn  =  n\  +  n2  +  wf  + — +n en  (7) 

Waarbij  het  superscript  e  de  eindtoestand  aanduidt.  Omdat  het  grootste  gedeelte 
van  een  eenheidscel  uit  vaste  stof  bestaat  wordt  ervan  uitgegaan  dat  in  de  eindtoe¬ 
stand  de  temperatuur  van  het  gas  gelijk  is  aan  de  temperatuur  van  de  gehele  een¬ 
heidscel  en  geldt  er: 

W+W+ . +  ^L+.nVL+W+ . +^L 

7j  R  T2R  TnR  TXR  T2R  TnR 


(8) 
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Na  vermeningvuldiging  met  R  en  buiten  haakjes  halen  van  Pe  wordt  de  volgende 
formule  voor  de  einddruk  in  de  buis  gevonden: 


r,v,  |  Mi  i. 


.+ 


PnV 2 


(9) 
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3  Computermodellen  en  werkzaamheden 


3.1  De  numerieke  code  ABAQUS 

Voor  de  berekeningen  is  de  eindige  elementen  code  ABAQUS  gebruikt.  Met  dit 
programma  kunnen  naast  sterkteberekeningen  ook  thermische  berekeningen  uitge- 
voerd  worden.  Naast  de  verschillende  type  randvoorwaarden  kunnen  interne 
warmtebronnen,  via  subroutines,  gedefinieerd  worden.  Voor  de  decompositie- 
warmte  van  een  explosieve  stof  is  in  het  verleden  reeds  een  speciale  ‘Fortran  77’- 
subroutine  geschreven.  Voor  een  uitvoerige  beschrijving  hiervan  wordt  naar  het 
rapport  ‘Thermische  berekeningen  aan  de  TNO-PML  Cook-off  test  in  ABAQUS 
met  de  Tarver-McGuire  reactiekinetiek’  [5].  De  decompositiereacties  hebben  een 
Arrhenius-structuur  waarbij  de  reactiesnelheid  afhankelijk  is  van  de  concentraties 
van  de  beginproducten  en  de  temperatuur.  Voor  de  eerste  modelverificatiebereke- 
ningen  met  HMX-PBX  zal  gebruikgemaakt  worden  van  de  zogenaamde  Tarver- 
McGuire-reactiekinetiek.  Hierbij  verloopt  de  decompositiereactie  via  meerdere 
endo-  en  exotherme  reactiesstappen.  Voor  deze  reactiekinetiek  zijn  de  waarden  uit 
de  literatuur  gebruikt.  De  decompositiereacties  van  de  AP/HTPB  en  AP/AN/PPG 
stuwstoffen  zijn  sterk  afhankelijk  van  allerlei  additieven  zodat  waarden  uit  de 
literatuur  niet  gebruikt  kunnen  worden.  Daarom  zijn  de  parameters  voor  deze 
stuwstoffen  via  thermogravimetriemetingen  (TGA)  bepaald. 


3.2  De  implementatie  van  het  drukmodel  in  het  Cook-off  model 
van  ABAQUS 

Als  uitgangspunt  voor  het  drukmodel  is  het  TNO-PML  Cook-off  model  genomen 
met  de  subroutine  HETVAL,  zoals  beschreven  staat  in  het  rapport  ‘Thermische 
berekeningen  aan  de  TNO-PML  Cook-off  test  in  ABAQUS  met  de  Tarver- 
McGuire  reactiekinetiek’  [5].  Naast  wat  kleine  veranderingen  in  de  invoerfile  voor 
de  rekencode  ABAQUS,  is  de  implementatie  van  het  drukmodel  hoofdzakelijk  te 
vinden  in  de  subroutine  HETVAL  (zie  bijlage  A).  In  het  uitgangsmodel  vindt  in 
deze  subroutine  de  berekening  van  de  reactiesnelheidsvergelijkingen  per  eenheids- 
cel  plaats,  gevolgd  door  een  berekening  van  de  warmteproductie/-verlies  ten  ge- 
volge  van  de  ontledingsreacties.  Met  de  berekende  gasproductieffactie  Fe  kan 
vervolgens  de  lokale  druk  in  de  eenheidscel  berekend  worden. 

Voordat  de  globale  druk  berekend  kan  worden  moet  nog  een  aantal  parameters 
bepaald  worden.  Omdat  de  subroutine  HETVAL  per  eenheids(reken)cel  aangeroe- 
pen  wordt  en  omdat  er  voor  de  berekening  van  de  globale  druk  Pe  een  som  van  een 
aantal  parameters  van  alle  eenheidscellen  berekend  moet  worden,  is  het  van  belang 
een  aantal  parameterwaarden  in  de  subroutine  te  kennen.  Voor  het  doorgeven  van 
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het  nummer  van  de  iteratieslag  van  de  berekening,  het  nummer  van  de  eenheidscel 

en  het  volume  van  deze  cel  is  een  aantal  regels  in  de  subroutine  toegevoegd. 

De  volgorde  van  berekenen  is  als  volgt. 

•  Nadat  de  gasproductiefactor  Fe  bepaald  is  (met  behulp  van  de  aanroep  van  de 
subroutine  odeint(..),  zie  bijlage  A  onder  het  kopje  ‘aanroep  van  het  oplossen 
van  de  diff  vlg....’),  worden  de  waarden  van  de  iteratieslag,  het  nummer  en  het 
volume  van  de  eenheidscel  opgehaald  (kopjes  ‘berekeningen  voor  node  en  ele- 
menten’  en  ‘ophalen  van  volume  van  eenheidscel’). 

•  Daama  wordt  de  lokale  druk  van  de  betreffende  eenheidscel  berekend  via  een 
subroutine  KPRESS. 

•  Daama  wordt  gecheckt  of  de  volgende  iteratieslag  van  de  berekening  al  gestart 
is.  Is  dit  wel  het  geval  dan  wordt  de  waarde  van  de  globale  druk  weggeschreven 
in  een  file.  Indien  de  waarde  van  de  globale  druk  groter  is  dan  400  MPa,  zal  de 
druk,  met  een  aantal  andere  parameters,  als  functie  van  de  tijd  worden  wegge¬ 
schreven. 

•  Dit  wordt  gevolgd  door  de  berekening  van  de  som  van  de  parameters  die  nodig 
zijn  voor  de  globale  druktoename  (kopje  ‘ophogen  van  de  som  voor  de  globale 
drukberekening’).  De  subroutine  eindigt  bij  labelnummer  200  en  spring  terug 
naar  het  hoofdprogramma. 

In  tabel  1  is  een  opsomming  gegeven  van  de  gebruikte  parameters  in  het  compu- 

termodel  voor  HMX-PBX. 
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Tabel  1 :  Gegevens  voor  HMX-PBX  voor  computerberekeningen. 


Tarver-McGuire  Kinetiek: 

HMX-PBX  (85%  HMX) 

Decomp,  warmte  q,  (J/kg)  10"6  (factor  10  kleiner) 

-0,0418 

Decomp,  warmte  q2  (J/kg)  10-6 

1,254 

Decomp,  warmte  q3  (J/kg)  10"6 

5,016 

Freq.  factor  ln(Z1)  (1/s) 

48,7 

Freq.  factor  ln(Z2)(1/s) 

37,5 

Freq.  factor  ln(Z3)  (1/s) 

28,1 

Activeringsenergie/  gasconst.  E^R  (K)  10'3 

26,5 

Activeringsenergie/  gasconst.  E2/R  (K)  10“3 

22,2 

Activeringsenergie/  gasconst.  E^R  (K)  10"3 

17,2 

Soortelijke  warmte  Cp(J/kg-K) 

957,4+2.281 82*T 

Soortelijke  geleidingscoef.  X  (J/m-s-K) 

0,317 

Max.  theoretische  dichtheid  pth 

1,640 

Praktische  dichtheid  p 

1,6166±0,01 

Smeltpunt  (K) 

520 

Smeltwarmte  (J/kg)  10'5 

2,09 

BKW-parameters  en  drukvergeiijkingen 

BKW-ct 

0,5 

BKW-p 

0,284 

BKW-K 

11,85 

BKW-0 

10200 

Eindtoestand  vast-vloeistof  fractie  Xc 

0,32 

Moleculaire  gewicht  Mw  (mol) 

33,3 

Covolume  Eyk 

471 

Initiele  dichtheid  gas  pa°(kg/m3)  103 

1,29 

3.3  De  gridmodellen 

Voor  de  TNO-PML  Cook-off  buis  wordt  gebruikgemaakt  van  een  reeds  bestaand 
grid  zoals  in  figuur  B.l  (bijlage  B)  is  weergegeven.  Het  betreft  hier  een  langsdoor- 
snede  van  een  buis  van  50  cm  met  2  doppen  waarbij  in  het  rekenmodel  ten  gevolge 
van  symmetrieassen  slechts  een  vierde  deel  van  de  doorsnede  doorgerekend  hoeft 
te  worden. 

In  figuur  B.2  in  bijlage  B  is  het  gebruikte  model  van  de  mogelijke  raketmotor 
gegeven.  Het  gebruikte  raketmodel  heeft  een  oorsponkelijke  lengte  van  1550  mm 
waarvan  slecht  604  mm  gemodelleerd  is.  De  totale  straal  van  de  raket  is  127  mm. 
De  dikte  van  de  stalen  casing  bedraagt  2,2  mm,  van  de  ‘inhabitor’  2  mm  en  van  de 
liner  bedraagt  de  dikte  1  mm.  De  straal  van  de  binnenboring  is  30  mm.  Hierdoor  is 
de  dikte  van  de  stuwstof  in  de  hoofdmotor  91,8  mm.  In  het  model  heeft  de  ontste- 
ker  een  totale  lengte  van  166  mm  en  bedraagt  de  buitendiameter  40  mm.  Deze 
bezit  een  casing  van  1 ,35  mm  dikte  een  interne  thermische  protectie  van  HTPB  van 
2  mm  en  een  webdikte  van  de  propellant  van  3  mm.  De  ontsteker  is  via  een  ‘Save 
and  Arm  Device’  (SAD)  aan  de  raket  verbonden  (in  het  model  rood  aangegeven). 
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Deze  vormt  de  ‘bottle  neck’  voor  het  warmtetransport  in  de  richting  van  de  ontste- 
ker.  In  de  berekeningen  zal  de  warmtegeleidingscoefficient  van  dit  materiaal 
gevarieerd  worden  om  de  invloed  hiervan  op  het  opwarmprofiel  te  kunnen  bepalen. 

Het  derde  gridmodel  is  een  model  van  de  ‘igniter’  zoals  deze  gebruikt  is  een  slow 
Cook-off  test.  Het  tweedimensionale  model  bestaat  uit  een  deel  van  de  igniter.  De 
straal  van  de  stalen  mantel,  met  een  dikte  van  1,35  mm,  is  14,2  mm.  Hierin  bevindt 
zich  een  stuwstoflaag  met  een  dikte  van  3  mm  zodat  een  cilindervormige  ruimte 
met  een  straal  van  9,85  mm  overblijft,  waarin  zich  lucht  bevindt. 


3.4  Typeberekeningen 

Met  de  bovenstaande  beschreven  gridmodellen  zijn,  zoals  reeds  vermeld,  verschil- 
lende  typeberekeningen  uitgevoerd.  Allereerst  zijn  berekeningen  uitgevoerd  om 
het  nieuw  ontwikkelde  thermisch/chemisch  drukmodel  te  verifieren.  Hiema  zijn 
temperatuurprofielberekeningen  van  het  ontwikkelde  raketmodel  uitgevoerd.  De 
eerste  serie  berekeningen  zijn  uitgevoerd  onder  het  ‘Feasibility  project’,  de  eerste 
fase.  In  de  tweede  fase  zijn  karakteriseringsberekeningen  uitgevoerd  gevolgd  door 
herhalingsberekeningen  van  de  TNO-PML  Cook-off  buis  maar  met  een  lagere 
opwarmsnelheid.  Als  laatste  zijn  aan  de  hand  van  de  gegevens  van  stuwstoffen  uit 
fase  een,  berekeningen  uitgevoerd  met  het  slow  Cook-off  gridmodel.  In  de  volgen- 
de  paragrafen  worden  deze  berekeningen  uitvoerig  toegelicht. 

3.4.1  Verificatieberekeningen  aan  de  hand  van  HMX-PBX  model;  eerste 
fase 

Uit  eigen  ervaring  met  experimenten  blijkt  dat  in  vele  gevallen  in  een  Cook-off 
reactie  een  ‘responsie’  optreedt  voordat  zelfopwarming  in  de  explosieve  stof  wordt 
waargenomen  [4],  De  buis  barst  open  ten  gevolge  van  de  drukopbouw  die  tot- 
standkomt  door  de  omzetting  van  de  explosieve  stof  in  gas.  De  hoeveelheid  vaste 
stof  die  dan  omgezet  is  is  in  de  orde  grootte  van  1-5%  van  de  vaste  stof  [6].  Na  het 
experiment  worden  daarom  grote  hoeveelheden  explosief  materiaal  teruggevonden. 
Het  berekenen  van  de  drukopbouw  in  de  buis  kan,  naast  het  berekenen  van  de 
temperatuurprofielen,  uitwijzen  of  het  einde  van  de  Cook-off  reactie  het  gevolg 
van  een  drukopbouw  of  van  een  samengaande  runaway  met  een  drukopbouwreac- 
tie  is. 

In  deze  berekeningen  wordt  het  nieuw  ontwikkelde  thermisch/chemisch  drukmodel 
in  ABAQUS  aan  de  hand  van  HMX-PBX  berekeningen  geverifieerd.  Indien  bij  de 
computerberekening  de  druk  in  de  buis  de  breekdruk  van  de  testbuis  bereikt 
(240  MPa)  voordat  een  wegloopreactie  optreedt,  zal  het  einde  van  de  gesimuleerde 
test  het  gevolg  zijn  van  een  drukopbouw.  In  dit  geval  is  de  drukopbouw  het  bepa- 
lende  stopcriterium  voor  de  berekening.  Indien  de  breekdruk  optreedt  terwijl  de 
runaway  al  optreedt  zal  een  complex  systeem  ontstaan  waarbij  de  drukopbouw  en 
de  drukaflaat  ten  gevolge  van  scheurvorming  met  elkaar  zullen  wedijveren.  Indien 
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de  drukaflaat  groot  genoeg  is,  zal  de  reactie  eindigen  met  een  milde  explosie. 
Indien  de  drukopbouw  echter  sterker  groeit  dan  de  drukaflaat  bestaat  de  kans  op 
een  deflagratie-detonatie  overgang.  Deze  berekeningen  worden  vervolgens  met 
experimenten  vergeleken. 

3.4.2  Thermische  berekeningen  zonder  interne  bronterm  met  het  raket- 
motormodel;  eerste  fase 

Om  een  beter  inzicht  te  krijgen  in  de  temperatuurverschillen  zijn  voor  het  raket- 
model  zuiver  thermische  geleidingsberekeningen  uitgevoerd.  In  deze  berekeningen 
bedraagt  de  gesimuleerde  opwarmsnelheid  3,3  °C/uur  en  3  °C/minuut.  Verder  is  de 
geleidingscoefficient  van  het  materiaal  dat  de  ‘bottle-neck’  vormt  voor  de  warm- 
testroom,  in  de  richting  van  de  ontsteker  (de  SAD-cilinder),  gevarieerd.  Voor  de 
geleidingscoefficient  van  dit  metaal  zijn  respectievelijk  de  waarden  10,  45  en 
390  J/smK  gekozen.  De  middelste  waarde  komt  overeen  met  de  geleidbaarheid  van 
een  standaardstaalsoort,  de  laatste  is  de  geleidingscoefficient  van  koper,  een  goede 
geleider,  en  de  eerste  is  voor  het  geval  er  slechte  contactgeleiding  optreedt  tussen 
de  schoefdraden  van  de  verschillende  metalen.  In  deze  berekeningen  is  voor  de 
overige  contactgeleidingen  van  een  ideale  geleiding  uitgegaan.  Wei  zijn  in  het 
model  de  verschillende  liners  tussen  de  metalen  omhulsels  en  de  stuwstoffen,  die 
de  warmtestromen  in  de  richting  van  de  stuwstoffen  beperken,  meegenomen. 

3.4.3  TNO-PML  Cook-off  berekeningen;  tweede  fase 

Als  tussenfase  zijn  metingen  verricht  met  een  TNO-PML  Cook-off  buis,  gevuld 
met  een  kaliumzout,  om  de  buis  beter  te  karakteriseren.  In  combinatie  met  compu- 
terberekeningen  is  een  aantal  waarden  voor  de  verschillende  warmte- 
overdrachtscoefficienten  tussen  de  verschillende  materialen  bepaald.  Tevens  is  nog 
een  aantal  berekeningen  uitgevoerd  met  HMX-PBX’en  met  een  opwarmsnelheid 
van  0,01  °C/s  en  vergeleken  met  de  experimentele  waarden  uit  rapport  [10]. 

Met  de  reactiekinetiek  van  de  HTPB/AP  en  PPG/AP/AN-stuwstuffen  die  bepaald 
zijn  in  de  eerste  fase,  zijn  in  de  tweede  fase  berekeningen  met  het  TNO-PML 
Cook-off  buismodel  en  een  slow  Cook-off  testmodel  uitgevoerd.  Bij  deze  bereke¬ 
ningen  is  gebruikgemaakt  van  de  geleidingscoefficienten  die  uit  de  Cook-off 
experimenten  afgeleid  zijn.  De  berekeningsresultaten  zijn  vervolgens  vergeleken 
met  experimentele  waarden  uit  de  eerste  fase. 
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4  Resultaten  van  fase  1 


4.1  Verificatieberekeningen  met  HMX-PBX  in  TNO-PML  Cook¬ 
off  buis 

Met  het  computerprogramma  zoals  beschreven  in  hoofdstuk  3,  zijn  berekeningen 
uitgevoerd  en  deze  zijn  vergeleken  met  experimentele  resultaten.  Om  een  idee  te 
geven  van  de  temperatuurgradienten  in  de  TNO-PML  Cook-off  buis  vlak  voor  de 
thermische  wegloopreactie,  is  in  figuur  B.2  in  bijlage  B  een  temperatuurcontour- 
plot  van  een  simulatie  met  HMX-PBX  te  zien.  In  deze  figuur  is  te  zien  dat  de  de 
wegloopreactie  in  de  explosieve  stof  vlak  bij  de  wand  plaatsvindt. 

In  figuur  1  is  een  temperatuur-tijd  plot  van  een  experiment  met  HMX-PBX  gege- 
ven.  De  temperaturen  zijn  gemeten  in  het  centrum  van  de  buis  in  radiele  richting 
op  respectievelijk  0,0, 0,44,  0,87,  1,31  en  1,75  cm  (0,0=Tcentrum,  1,75  =  Twand) 
van  het  centrum  van  de  buis  (zie  ook  figuur  4).  Daamaast  is  nog  een  thermokoppel 
Tdop  in  de  dop  geplaatst. 


Figuur  1:  Temperaruur-tijd  curven  van  een  experiment  met  HMX-PBX  in  de  TNO-PML 

Cook-ojf  buis. 

De  opwarmsnelheid  in  dit  experiment  was  0,05  °CI s.  Na  een  zeer  kleine  afwijking 
van  de  lineaire  stijging  (vanaf  4000  s)  van  de  temperaturen  in  het  inwendige  van 
de  buis,  barst  de  buis  open  ten  gevolge  van  een  drukopbouw  van  de  gevormde 
gassen  bij  een  wandtemperatuur  van  229  °C.  De  afwijking  is  het  gevolg  van  endo- 
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therme  reacties  in  het  begin  van  de  ontledingsreactie  die  meestal  energie  kosten 
omdat  grote  ketens  opgebroken  moeten  worden. 

In  figuur  2  zijn  de  resultaten  van  de  temperatuurberekeningen  aan  deze  buis  gege- 
ven.  In  figuur  3  is  de  bijbehorende  globale  drukopbouw  weergegeven. 


Tijd  (s) 


Figuur  2:  Resultaten  van  de  berekening  van  de  temperatuur-tijd  curven  van  de  TNO- 

PML  Cook-off  buis  gevuld  met  HMX-PBX. 
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Tijd  (s) 


Figuur  3:  Druk-temperatuur-tijd  curven  van  TNO-PML  Cook-off  buis  met  HMX-PBX. 

Volgens  de  simulatie  vindt  de  thermische  wegloopreactie  plaats  bij  een  wandtem- 
peratuur  van  231  °C.  De  druk  bereikt  in  de  simulatie  na  4366  s  de  statische  breek- 
druk  van  240  MPa  van  de  buis  bij  een  temperatuur  van  227  °C.  Dit  is  4  graden 
voordat  de  thermische  wegloopreactie  ingezet  wordt.  Aangezien  in  het  experiment 
de  buis  ten  gevolge  van  een  drukopbouw  gebarsten  is  en  nog  grote  resten  niet- 
gereageerd  HMX-PBX  teruggevonden  zijn,  blijkt  de  berekening  goed  met  de 
werkelijkheid  overeen  te  komen.  Hieruit  volgt  dat  het  model  een  goede  omschrij- 
ving  van  de  werkelijkheid  kan  geven;  namelijk  een  berekende  breekdruk  met  een 
corresponderende  wandtemperatuur  van  227  ±  2  °C  ten  opzichte  van  een  gemeten 
wandtemperatuur  van  229  ±  2  °C  waarbij  de  buis  bezwijkt  ten  gevolge  van  een 
drukopbouw  van  de  gevormde  gassen. 


4.2  Thermische  geleidingsprofielen  van  de  raketmotor 

Na  de  Cook-off  HMX-PBX  berekeningen  zijn  de  temperatuurcontourprofielen  van 
de  raketmotor  bij  twee  opwarmsnelheden  en  drie  verschillende  geleidingscoeffi- 
cienten  van  de  SAD-cilinder  berekend.  In  bijlage  B  zijn  de  resultaten  van  deze 
berekeningen  gegeven.  De  eerste  drie  temperatuurprofielen  1 1,  12  en  13  zijn  de 
resultaten  van  de  berekeningen  met  een  opwarmsnelheid  van  0,05  °C/s.  De  contou- 
ren  geven  de  toestand  aan  om  het  moment  dat  de  casing  een  temperatuur  heeft  van 
300  °C.  Het  profiel  van  de  opwarming  is  dan  nog  niet  constant  maar  volgens  de 
experimenten  [1]  is  er  in  het  geval  met  een  explosieve  stof  reeds  een  thermische 
wegloopreactie  opgetreden. 
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Indien  we  naar  de  drie  opeenvolgende  contouren  kijken  is  zichtbaar  dat  bij  de 
laagste  waarde  van  de  geleidmgscoefficient  de  warmte  het  minst  diep  in  de  ontste- 
ker  dringt.  Deze  heeft  nauwelijks  een  temperatuur  van  50  °C  terwijl  de  tempera- 
tuur  in  de  wand  en  in  hoofdstuwstof  al  een  temperatuur  van  300  °C  heeft  bereikt. 
Bij  de  twee  opeenvolgende  geleidingscoefficienten  van  45  en  390  J/ms2  K  heeft  de 
ontsteker  al  een  temperatuur  van  meer  dan  100,  respectievelijk  150  °C  bereikt. 
Toch  zal  bij  deze  opwarmsnelheid  de  hoofdlading  eerder  ontsteken  dan  de  ontste¬ 
ker  daar  de  verschillen  te  groot  zijn  om  de  ontsteker  eerder  te  laten  ontsteken. 
Omdat  er  slechts  een  fractie  van  de  hoofdlading  op  hoge  temperatuur  is  en  het 
overgrote  deel  zich  nog  op  een  relatief  lage  temperatuur  bevindt,  zullen  de  me- 
chanische  eigenschappen  hoogstwaarschijnlijk  nog  zo  goed  zijn  dat  een  redelijk 
normale  verbranding  van  de  raket  verwacht  kan  worden.  In  figuur  B.10  in 
bijlage  B  is  aan  de  hand  van  een  indicatieve  berekening  te  constateren  dat  de 
wegloopreactie  start  in  de  rechterbovenhoek  van  het  raketmodel. 

In  de  drie  daaropvolgende  contouren  zijn  de  eindprofielen  van  de  simulatie  weer- 
gegeven  bij  een  opwarmsnelheid  van  3,3  °C/uur.  De  contouren  zijn  de  tempera- 
tuurprofielen  op  de  tijdstippen  dat  de  casing  een  temperatuur  van  200  °C  heeft.  In 
deze  toestand  heeft  zich  een  constant  temperatuurprofiel  gevormd  zodat  bij  verdere 
opwarming  geen  veranderende  verschillen  in  temperatuur  meer  zullen  optreden 
totdat  er  zelfopwarming  van  de  explosieve  stof  gaat  optreden. 

In  alle  contouren  bij  een  opwarmsnelheid  van  3,3  °C/uur  bedraagt  het  maximale 
temperatuurverschil  12  °C.  Bij  de  verschillende  geleidingscoefficienten  van  10,  45 
en  390  J/ms2  K  van  de  S  AD-cilinder  zijn  de  temperaturen  in  de  ontsteker  respec¬ 
tievelijk  188,  191  en  196  °C,  terwijl  de  temperatuur  van  de  casing  200  °C  be¬ 
draagt. 
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5  Resultaten  van  fase  2 


5.1  Inleiding 

In  het  kader  van  het  project  ‘Thermische  Initiatie’  is  een  aantal  simulaties  in 
combinatie  met  experimenten  uitgevoerd.  Allereerst  is  een  aantal  zogenaamde 
karakteriseringsexperimenten  uitgevoerd  om  de  warmteoverdrachtscoefficienten 
van  verschillende  overdrachten  in  de  TNO-PML  Cook-off  buis  te  bepalen.  Tevens 
zijn  HMX-PBX  simulaties  van  de  TNO-PML  Cook-off  buis  uitgevoerd  maar  nu 
met  een  opwarmsnelheid  van  0,01  °C/s,  hetgeen  een  factor  5  langzamer  is  dan  de 
simulaties/experimenten  uit  fase  1.  Bij  deze  opwarmsnelheid  zal  de  thermische 
wegloopreactie  vanuit  het  midden  plaatsvinden.  Deze  simulaties  zijn  wederom 
vergeleken  met  de  experimentwaarden.  Als  laatste  zijn  simulaties  uitgevoer  met 
een  slow  Cook-off  model  en  vergeleken  met  de  experimentele  waarden.  De  ge- 
bruikte  gegevens  van  zijn  afkomstig  van  de  ‘feasibility  study  of  an  IM  ignitor’  [3]. 
Het  betreft  hier  nulde  orde  en  niet-nulde  orde  reactiekinetiekparameters. 


5.2  Karakteriseringsexperimenten  versus  berekeningen 

Het  doel  van  de  karakteriseringsexperimenten  was  de  TNO-PML  Cook-off  buis 
thermisch  te  karakteriseren,  hetgeen  inhoudt  dat  onder  meer  de  randvoorwaarden 
en  de  warmteoverdrachtscoefficienten  bepaald  worden.  Hiervoor  zijn  op  verschil¬ 
lende  overgangen  van  het  ene  naar  het  andere  materiaal  thermokoppels  geplaats  en 
is  een  bepaald  temperatuurtraject  gevolgd.  Met  hetzelfde  traject  zijn  eveneens 
simulatie  uitgevoerd  en  zijn  de  verschillende  warmteoverdrachtscoefficienten 
gevarieerd.  Uit  de  vergelijking  van  de  resultaten  van  deze  metingen  met  de  simu¬ 
laties  is  een  schatting  voor  de  verschillende  warmteoverdrachtscoefficienten 
gemaakt. 
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Figuur  4:  Ingekorte  schematische  tekeningen  van  de  TNO-PML  Cook-off  buis  met 

locatie  van  de  thermokoppels  en  overdrachtsgebieden  ( dikke  lijnen)  in  de  ka- 
rakteriseringsmetingenf  -berekeningen. 
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In  figuur  4  is  een  schematische  tekening  van  de  TNO-PML  Cook-off  buis  gegeven 
met  de  locaties  van  de  thermokoppels.  Eveneens  zijn  in  deze  figuur  de  namen  van 
de  bepaalde  warmteoverdrachtscoefficienten  via  een  dikke  lijn  aangegeven.  Gap- 
leng  is  de  overdracht  tussen  de  teststof  en  de  stalen  buis,  gapschr  tussen  de 
schroefdraad  van  de  dop  en  de  buis,  gapcop  is  de  warmteoverdracht  tussen  de  dop 
en  de  buis  via  de  koperen  cilinder  en  gapsple  de  overdracht  tussen  de  dop  en  de 
teststof.  Daamaast  zijn  er  nog  warmteoverdrachtscoefficienten  aan  de  lucht  via  de 
cilindrische  kant  van  de  dop  ( dopcyl)  en  de  vlakke  kant  van  de  dop  (dopvlak).  In 
figuur  C.l  in  bijlage  C  is  een  resultaat  van  een  karakteriseringsmeting  gegeven.  De 
nummers  van  de  thermokoppei  komen  hierbij  overeen  met  de  nummer  in  de  teke¬ 
ning.  In  figuur  C.2  in  bijlage  C  is  een  simulatie  van  dit  experiment  gegeven  nadat 
de  warmteoverdrachtscoefficienten  bepaald  zijn.  De  volgende  waarden  zijn  gevon- 
den  (tabel  2). 

Tabel  2:  Warmteoverdrachtscoefficienten  van  de  TNO-PML  Cook-off  buis. 


(J/s-m2-K) 

Niet-ge'fsoleerd 

Gei'soleerd  (6  cm) 

Gapleng 

350 

350 

Gapschr 

1000 

1000 

Gapcop 

1200 

1200 

Gapsple* 

350,  6,0 

350,  6,0 

dopcyl 

5,0,  6,0  (afhank.  van  T) 

0,5 

dopvlak 

6,0,  7,7  (afhank.  van  T) 

0,5 

*  Indien  er  direct  contact  met  de  dop  is  geldt  de  hoge  waarde,  is  er  een  luchtspleet  tussen  de  teststof 
en  de  dop  dan  geldt  de  lage  waarde. 


5.3  HMX-PBX  berekeningen  versus  experiment 

Van  de  gegevens  uit  enerzijds  de  karakteriseringsexperimenten/-berekeningen  en 
anderzijds  uit  de  berekeningen  van  fase  1,  is  gebruikgemaakt  bij  de  berekening  van 
de  simulaties  met  HMX-PBX  met  een  opwarmsnelheid  van  0,01  °C/s.  Deze  bere¬ 
keningen  zijn  wederom  vergeleken  met  de  experimentele  waarden. 
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Figuur  5:  Resultaat  van  een  HMX-PBX  simulatie  met  een  opwarmsnelheid  van  0,01  °C/s. 

In  figuur  5  is  het  resultaat  van  een  HMX-PBX  simulatie  gegeven  met  daamaast  de 
uitvergroting  van  de  laatste  fase.  De  verticale  donkere  lijn  geeft  hierbij  het  moment 
aan  waarbij  de  breekdruk  van  de  Cook-off  buis  overschreden  wordt  en  daarmee  het 
einde  van  de  simulatie.  In  figuur  6  is  een  experimenteel  resultaat  met  dezelfde 
opwarmsnelheid  gegeven.  Tevens  is  ook  van  dit  experiment  een  vergroting  van  de 
eindfase  van  het  experiment  gegeven.  Een  vergelijking  van  de  vier  grafieken  laat 
duidelijk  zien  dat  de  berekeningen  goed  overeenkomen  met  het  experiment  tot  en 
met  de  eindfase  toe.  Naast  de  temperatuurverschillen  in  de  buis  is  een  wegloopre- 
actie  ten  gevolge  van  de  ontleding  van  de  explosieve  stof  in  beide  figuren  waar  te 
nemen. 


Tlme(s) 


Figuur  6:  Het  resultaat  van  een  experiment  met  HMX-PBX  in  de  TNO-PML  Cook-off  buis  met  een  opwarmsnel¬ 

heid  van  0,01  °Cls. 
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5.4  Slow  Cook-off  berekeningen  versus  het  experiment 

In  de  laatste  fase  zijn  met  de  gegevens  van  de  ‘Feasibility  study  of  an  IM-Ignitor’, 
berekeningen  uitgevoerd  met  een  slow  Cook-off  model  zoals  beschreven 
hoofdstuk  3.  In  de  evaluatie  van  het  ‘feasibility-project’  [3]  zijn  de  gegevens 
samengevat. 

In  tabellen  3,  4, 5  en  6  is  een  overzicht  gegeven  van  de  reactiekinetiekparameters 
die  van  belang  zijn  voor  de  berekeningen.  De  waarden  zijn  afkomstig  van  thermo- 
gravimetrieanalyses  (TGA)  en  DSC-metingen. 


Tab  el  3:  Resultaten  van  de  kinetiekberekeningen  voor  een  AP/HTPB-stuwstof 


Snelheid 

Orde 

-E/r 

Ea  (kJ/mol) 

In  (A) 

A  (min/1) 

1  °C/min. 

1 

-6,75975 

56,20 

10,541 

3,78E4 

2 

-8,60023 

71,21 

14,665 

2,34E6 

0,5 

-5,83951 

48,55 

8,4792 

4,81  E3 

Tabel  4: 

Resultaten  voor  de  reactiekinetiek  voor  een  APlHTPB-stuwstof  (nulde  orde). 

Conversie  2% 

Ea=1 87,42  kJ/mol 

Qo=5,4515E19  J/kg-s 

Conversie  4% 

Ea=232,04  kJ/mol 

Qo=4,2952E23  J/kg-s 

Conversie  8% 

Ea=118,22  kJ/mol 

Qo=1,1543E11  J/kg-s 

Tabel  5: 

Resultaten  van  de  kinetiekberekeningen  voor  een  PPG/AP/AN-stuwstof 

Snelheid 

Orde 

-E/r 

Ea  (kJ/mol) 

In  (A) 

A  (min/1) 

1  °C/min. 

1 

-6,81471 

56,66 

8,757 

6,36E3 

2 

-8,29090 

68,93 

12,210 

2,01  E5 

0,5 

-6,07666 

50,52 

7,0312 

1,1 3E3 

Tabel  6: 

Resultaten  voor  de  reactiekinetiek  voor  een  PPG/AP/AN-stuwstof  (nulde 

orde). 

Conversie  1% 

Ea=1 07,64  kJ/mol 

Qo=5,51  El  1  J/kg-s 

Conversie  2% 

Ea=1 30,42  kJ/mol 

Qo=8,45  El 3  J/kg-s 

Conversie  4% 

Ea=1 80,50  kJ/mol 

Qo=1,47  El  9  J/kg-s 

Omdat  de  frequentiefactoren  (ln(A»  voor  beide  stuwstoffen  vrij  laag  waren  is 
gekeken  of  de  berekening  van  de  parameters  ook  overeenkwamen  met  de  gemeten 
TGA-waarden.  In  figuur  7  is  het  resultaat  van  een  ontledingsberekening  voor  een 
PPG/AP/AN-stuwstof  gegeven  samen  met  de  origineel  gemeten  TGA-waarden. 

Tot  een  omzetting  van  20%  (l-a=0,8  in  de  grafiek)  komt  de  gemeten  waarde  goed 
overeen  met  de  berekende  waarde.  Omdat  bij  omzettingen  van  minder  dan  5% 
explosieve  stof  een  ‘runaway’-reactie  al  op  voile  gang  is,  zou  dit  gebied  voldoende 
moeten  zijn  voor  een  simulatie  van  een  wegloopreactie  in  een  Cook-off  bereke¬ 
ning. 
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Figuur  7:  Omzetting  van  een  AP/ANIPPG  stuwstof  in  een  TGA  met  een  opwarmsnel - 

heid  van  1  °C/min .  samen  met  de  berekende  curve  voor  de  Cook-off  simula- 
tie . 

In  tabel  7  is  een  overzicht  gegeven  van  de  experimentele  en  berekende  waarde  van 
de  Cook-off  experimenten  met  stuwstoffen.  Zoals  uit  de  tabel  blijkt,  zijn  er  geen 
waarden  voor  de  berekeningen  met  niet-nulde  orde  reacties.  In  de  berekeningen 
start  de  decompositie  van  de  stuwstof  al  rond  de  150  °C  en  heeft  een  volledige 
omzetting  van  de  stof  al  bij  temperaturen  rond  de  240  °C  plaatsgevonden.  Even- 
eens  was  de  waarde  voor  de  totale  energie-inhoud  klaarblijkelijk  zo  laag  dat  geen 
‘runaway*  in  de  berekening  geconstateerd  werd. 

De  nulde  orde  kinetiek  doet  het  daarentegen  beter  wat  betreft  de  AP/HTPB- 
stuwstof.  Bij  een  omzetgraad  van  2%  is  de  wandtemperatuur  op  het  moment  van 
de  thermische  wegloopreactie  Twand=267  °C,  en  is  de  afwijking  ten  opzichte  van 
het  experiment  (272  °C)  5  °C.  Bij  een  omzetting  van  4%  is  de  afwijking  5  °C  naar 
boven  en  is  de  berekende  explosietemperatuur  277  °C.  Voor  de  AP/ AN/PPG 
stuwstof  komen  de  waarden  niet  overeen  met  de  experimentele  waarden. 


Tabel  7:  Overzicht  van  de  resultaten  van  de  slow  Cook-off  simulate  en  experimenten . 


Stuwstof 

Texpl  experimenteel 

Texpl  Simulatie 
nulde  orde 

Texpl  Simulatie 
A->B  orde  n=0,5 

AP/HTPB 

272 

267  (2%) 

geen 

277  (4%) 

geen 

AP/AN/PPG 

179 

267  (2%) 

geen 

265  (4%) 

geen 
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6  Discussie 


Thermische  berekeningen  met  HMX-PBX 

Indien  de  zuivere  thermische  simulatieresultaten  van  HMX-PBX  met  de  experi- 
mentele  resultaten  vergeleken  worden,  kan  gesteld  worden  dat  de  licht-aangepaste 
Tarver  McGuire-kinetiek  met  een  driestaps  reactiekinetiekmodel  goed  voldoet.  Bij 
de  opwarmsnelheid  van  0,05  °C/s  en  van  0,01  °C/s  zijn  de  verschillen  voor  de 
explosietemperatuur,  de  wandtemperatuur  op  moment  van  de  ‘explosie’,  in  verge- 
lijking  tot  de  experimenten,  niet  groter  dan  enkele  graden.  Uit  de  vergelijking  van 
de  temperatuurgradienten  van  beide  experimenten  met  de  resultaten  van  de  simu- 
laties  kan  ook  gesteld  worden  dat  de  berekeningen  goed  met  de  praktijk  overeen- 
komen. 

Duidelijk  is  echter  dan  bij  het  experiment  met  een  opwarmsnelheid  van  0,05  °C 
nog  geen  duidelijke  wegloopreactie  te  zien  is,  maar  slechts  een  lichte  afwijking 
van  de  lineaire  opwarmsnelheid.  Tevens  kan  uit  de  fragmentatie  van  de  buis  afge- 
leid  worden  dat  er  slechts  sprake  is  van  het  openbreken  van  de  buis  ten  gevolge 
van  een  verhoogde  druk.  In  enkele  experimenten  met  PBX-en  zijn  zelfs  na  het 
experiment  nog  grote  delen  niet  gereageerde  explosieve  stof  temggevonden.  Bij 
het  langdurige  experiment  met  een  opwarmsnelheid  van  0,01  °C/s  is  er  wel  sprake 
van  het  begin  van  een  wegloopreactie,  vanuit  het  midden  van  de  buis,  maar  is  deze 
nog  niet  maximaal  qua  snelheid. 

Het  gecombineerde  thermisch/chemisch/drukmodel 

Door  gebruik  te  maken  van  het  drukmodel  kan  blijkbaar  een  betere  schatting  van 
het  moment  van  reactie  gemaakt  worden.  Bij  het  0,05  °C/s  experiment  blijkt  dui¬ 
delijk  uit  de  simulatie,  net  als  in  het  experiment,  dat  de  interne  druk  de  breekdruk 
van  de  buis  overschrijdt  voordat  de  thermische  wegloopreactie  plaatsvindt.  Ook  in 
de  simulatie  van  het  langzame  experiment  wordt  de  breekdruk  van  de  buis  bereikt 
voordat  de  wegloopreactie  maximaal  is. 

In  hoeverre  het  drukmodel  overeenkomt  met  de  praktijk  is  moeilijk  te  zeggen.  Uit 
onderzoek  [11]  blijkt  dat  het  model  niet  direct  uit  de  lucht  komt  vallen.  Uit  de 
praktijk  blijkt  dat  er  inderdaad  holten  gevormd  worden  waarin  de  decompositie- 
gassen  zich  verzamelen  waardoor  de  holten  groeien.  De  gebruikte  gasfase- 
vergelijking  is  slechts  een  model  van  een  toestand  zoals  deze  in  een  prefase  van  de 
Cook-off  reactie  verwacht  wordt.  De  gebruikte  constante  parameters  zullen  zeker 
niet  voor  het  gehele  Cook-off  proces  gelden.  Eveneens  zal  de  mate  van  poreusiteit 
en  de  verandering  hiervan  als  functie  van  de  ontledinggraad,  en  de  diffusiesnelheid 
van  het  gas,  een  belangrijke  rol  spelen  in  berekening  van  de  lokale  en  globale  druk. 
Ondanks  het  groot  aantal  randvoorwaarden  lijkt  het  er  op  dat  met  behulp  van  de 
drukberekeningen  een  beter  stopcriterium  is  gevonden  voor  de  numerieke  simula¬ 
tie  van  een  Cook-off  experiment.  Voor  de  bevestiging  zullen  nog  enkele  zeer 
‘goede’  drukmetingen  nodig  zijn. 
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Mede  dankzij  de  karakteriseringstesten/-simulaties  kunnen  goede  schattingen  van 
de  temperatuurprofielen  van  de  TNO-PML  Cook-off  buis  gemaakt  worden  met 
behulp  van  de  computercode.  Alle  berekende  temperaturen  verschillen  slechts 
enkele  graden  van  de  gemeten  waarden,  met  inbegrip  van  de  uiterste  doptempera- 
turen.  De  belangrijkste  warmteoverdrachtscoefficienten  zijn  bekend  en  zullen 
bruikbaar  zijn  in  toekomstige  berekeningen  met  de  TNO-PML  Cook-off  buis. 

De  berekeningen  van  temperatuurprofielen  van  het  raketsysteem 
De  raketsysteemtemperatuurprofielberekeningen  geven  een  goede  indruk  van  de 
temperatuurgradienten  die  verwacht  kunnen  worden  bij  dergelijke  systemen.  Het 
blijkt  dat  bij  een  intermediaire  opwarmsnelheid  van  3  °C/min.  enorme  tempera- 
tuurverschillen  gevonden  worden  tussen  het  interne  van  de  raket  en  de  buitenwand. 
Tevens  blijkt  dat  bepaalde  doorgangen  een  enorme  invloed  hebben  op  de  warmte- 
doorstroom  in  de  richting  van  de  ontsteker.  Er  werden  temperatuuurverschillen 
gevonden  van  250  tot  150  °C. 

Ook  bij  een  opwarmsnelheid  van  3,3  °C/uur  (slow  Cook-off )  zijn  nog  tempera- 
tuurverschillen  tot  12  °C  gevonden.  Hierbij  hebben  de  doorgangen  wel  invloed  op 
de  maximale  verschillen,  maar  zijn  de  verschillen  beduidend  minder.  De  verschil¬ 
len  bedroegen  4  tot  12  °C.  In  deze  analyses  is  opwarming  van  de  stuwstoffen  niet 
meegenomen. 

Om  de  analyse  meerwaarde  te  geven  is  gekeken  hoe  de  stoffen  reageren  in  de  slow 
Cook-off,  SCB  en  TNO-PML  Cook-off  testen  [3].  In  een  niet  opgesloten  toestand 
reageert  de  HTPB/AP-stuwstof  bij  ongeveer  272  °C  terwijl  de  PPG/AP/AN- 
stuwstof  bij  179  °C  een  reactie  geeft.  In  opgesloten  toestand,  zonder  vrije  ruimte  in 
het  testitem,  ontlopen  ze  elkaar  niet  veel.  Bij  een  opwarmsnelheid  van  0,05  °C/s 
geeft  de  HTPB/AP-stuwstof  bij  200  °C  een  reactie  terwijl  de  PPG/AP/AN-stuwstof 
bij  190  °C  een  reactie  geeft.  De  volgende  denkbare  cases  kunnen  zich  nu  voor- 
doen. 

Case  1 

Indien  de  raketmotor  als  een  gesloten  systeem  beschouwd  wordt,  terwijl  de  ontste¬ 
ker  zijn  gassen  kwijt  kan  in  de  centrale  boring  van  de  raketmotor,  kan  volgens  de 
experimentele  gegevens  het  volgende  van  toepassing  zijn. 
a  Indien  alle  stuwstof  van  HTPB/AP  gemaakt  is  zal  de  hoofdlading  vroeger  dan 
de  ontsteker  ontbranden  en  kan  ten  gevolge  van  de  slechte  mechanische  eigen- 
schappen  bij  deze  hoge  temperatuur  een  zeer  heftige  reactie  ontstaan. 
b  Indien  de  hoofdlading  van  HTPB/AP  en  de  ontsteker  van  PPG/AP/AN  ge¬ 
maakt  zijn,  zal  de  hoofdlading  bij  een  lagere  temperatuur  dan  200  °C  ontsteken 
(bij  lage  opwarmsnelheid  van  3,3  °C/uur  zal  in  opgesloten  toestand  de  weg- 
loopreactie  altijd  lager  zijn  dan  bij  een  opwarmsnelheid  van  3  °C/min.)  terwijl 
de  PPG  stuwstof  bij  179  °C  ontsteekt.  Indien  de  ontsteker  eerder  en  op  een 
normale  manier  brandt  dan  de  hoofdlading,  bepaalt  het  brandgedrag  van  de 
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hoofdlading  bij  deze  temperatuur  of  de  verbranding  overgaat  in  een  heftige  re- 
actie.  Indien  de  hoofdlading  eerder  ontbrandt  dan  de  ontsteker  zal  er  geen  ver- 
schil  zijn  met  case  la. 

Case  2 

Indien  de  raketmotor  ook  als  open  systeem  beschouwd  wordt  zal  het  volgende 
gelden. 

a  Indien  alle  stuwstof  van  HTPB/AP  gemaakt  is  zal  de  raket  vanuit  het  midden 
van  de  hooflading  tot  ontbranding  komen.  Net  als  in  case  la  zal  dit  hoogst- 
waarschijnlijk  tot  een  heftige  reactie  leiden  ten  gevolge  van  de  verslechterde 
mechanische  eigenschappen. 

b  Indien  de  ontsteker  van  PPG-stuwstof  gemaakt  is  zal  deze  hoogstwaarschijn- 
lijk  eerder  ontsteken  dan  de  HTPB-stuwstof  daar  de  ontbrandingstemperaturen 
van  de  beide  stuwstoffen  in  een  open  slow  Cook-off  test  sterk  uit  elkaar  liggen. 
Hierdoor  zal,  indien  de  ontsteker  normaal  functioneert  en  de  HTPB-stuwstof 
nog  een  normaal  brandgedrag  vertoont,  de  raketmotor  redelijk  normaal  kunnen 
functioneren. 

Case  3 

Indien  beide  systemen  als  gesloten  beschouwd  worden  kan  het  volgende  gelden. 
a  Indien  alle  stuwstof  van  HTPB/AP  gemaakt  is  zal  net  zoals  in  case  la  de 
raketmotor  hoogstwaarschijnlijk  heftig  reageren. 
b  Indien  de  ontsteker  van  PPG-stuwstof  gemaakt  is  zullen  beide  elkaar  niet  veel 
ontlopen  qua  ontbrandingstemperatuur.  Door  de  grote  bulk  stuwstof  in  de 
hoofdlading  is  het  waarschijnlijk  dat  de  gegenereerde  warmte  hier  niet  kan 
wegebben  waardoor  de  ontsteking  in  het  midden  van  de  hoofdlading  zal  be- 
ginnen  voordat  de  ontsteker  ontbrandt.  Hierdoor  zal  weer  een  reactie  optreden 
zoals  in  case  la. 

Door  het  feit  dat  niet  goed  bepaald  kan  worden  of  in  een  raketmotor  een  breek- 
membraan  en  de  boring  voor  een  gesloten  dan  wel  open  systeem  zorgen,  terwijl  de 
reacties  van  de  beide  stuwstoffen  in  een  open  en  gesloten  systeem  volledig  ver- 
schillend  reageren  kan  geen  eenduidig  antwoord  gegeven  worden  hoe  aan  de  hand 
van  de  temperatuurprofielen  de  raketmotor  bij  langzame  opwarming  zal  reageren. 

Slow  Cook-off  berekeningen 

Voor  de  slow  Cook-off  berekeningen  komen  de  berekende  reactiesnelheidsverge- 
lijkingen  niet  overeen  met  de  experimentele  waarden.  Wellicht  is,  door  beperkte 
hoeveelheid  tijd  die  hieraan  besteed  kon  worden,  niet  genoeg  onderzoek  verricht 
naar  het  bepalende  reactiemechanisme  en  overall  bepalende  reactiestap  bij  de 
ontleding  van  de  stuwstoffen.  Bij  de  ontleding  van  complexe  composities  zoals  een 
HTPB/AP  [12,  13,  14]  met  toevoegingen  die  de  verbranding  moeten  regelen,  zijn 
zoveel  concurrerende  reactiestappen  dat  een  diepgaander  onderzoek  nodig  is  om 
de  belangrijke  stappen  te  bepalen  en  te  kwantificeren.  Tevens  is  het  belangrijk  om 
te  weten  of  de  kinetische  parameters  met  een  open  dan  wel  gesloten  systeem 
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bepaald  zijn.  Bij  de  ontleding  van  vele  energetische  materialen,  zoals  AP,  komen 
gassen  vrij  die  een  katalytische  werking  kunnen  hebben  op  de  reactiesnelheid. 
Indien  deze  in  een  open  systeem  kunnen  vluchten  zal  dit  wel  degelijk  invloed 
hebben  op  de  bepaalde  kinetiek  die  van  toepassing  moet  zijn  voor  een  gesloten 
systeem. 

Een  andere  reden  voor  de  niet-correcte  bepaling  van  de  kinetiekparameters  is  de 
volgende.  Bij  vele  energetische  materialen  bestaat  de  ontleding  van  het  materiaal 
uit  een  aantal  al  dan  niet  parallele  reactiestappen,  waarbij  de  eerste  stap  een  endo- 
therme  stap  is,  waarbij  de  grote  ketens  in  kleine  deelproducten  opgesplitst  worden. 
Deze  stap  wordt  veelal  gevolgd  door  licht-exotherme  stappen  met  als  laatste  een  of 
meer  sterk-exotherme  stappen  die  verantwoordelijk  zijn  voor  de  grote  hoeveelhe- 
den  warmte  die  vrijkomen.  Bij  een  TGA  echter  wordt  niet  naar  de  energietoena- 
men  of  -afname  gekeken  maar  naar  de  afname  van  de  massafractie  als  functie  van 
de  temperatuur.  Aan  de  hand  van  deze  afname,  in  de  eerste  fase  van  de  ontleding, 
is  de  bepaling  van  de  kinetiek  uitgevoerd.  Klaarblijkelijk  is  dit  niet  de  belangrijk- 
ste  fase,  waarbij  de  grote  hoeveelheden  warmte  vrijkomen,  voor  de  ontleding  van 
deze  stuwstoffen.  Dus,  ondanks  dat  de  massafractieafname  goed  gemodelleerd 
wordt  in  de  beginfase  van  de  ontleding,  hoeft  dit  niet  te  betekenen  dat  dit  een 
goede  modellering  is  voor  de  ontleding  van  de  stuwstoffen. 
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7  Conclusies 


Er  is  een  aantal  computermodellen  ontwikkeld  waarmee  ofwel  eenvoudige  ther- 
misch-chemische  berekeningen  ofwel  meer  complexe  gecombineerde  thermisch- 
chemische  drukberekeningen  uitgevoerd  kunnen  worden.  Naast  gegevens  uit  de 
literatuur,  zijn  thermische-  en  reactiesnelheidsparameters  bepaald,  die  nodig  zijn 
voor  deze  modellen.  Met  de  ontwikkelde  computerprogramma’s  zijn  vervolgens 
TNO-PML  Cook-off  buis,  raketsysteem  en  slow  Cook-off  berekeningen  uitgevoerd 
en  vergeleken  met  de  experimentele  waarden. 

Met  de  resultaten  van  de  temperatuurprofielberekeningen  van  het  raketsysteem  is 
een  beter  inzicht  verkregen  in  het  warmtetransport  door  het  systeem.  Dit  kan  bij  de 
ontwikkeling  van  nieuwe  systemen  een  bijdrage  leveren  aan  het  ontwerp.  Een 
betere  inschatting  van  een  mogelijke  responsie  van  een  ratketsysteem  op  een 
thermische  stimulus  kan  gemaakt  worden  indien  naast  deze  berekeningen  de  resul¬ 
taten  van  de  verschillende  Cook-off  experimenten  toegevoegd  worden. 

Het  blijkt  dat  de  reactiesnelheidsvergelijkingen  die  bepaald  zijn  voor  de  stuwstof- 
fen  HTPB/AP  en  PPG/AP/AN  niet  overeenkomen  met  de  werkelijkheid.  Hiervoor 
is  een  aantal  redenen  genoemd  waaruit  blijkt  dat  de  bepaling  van  de  reactiesnel- 
heidsvergelijking  van  een  energetisch  materiaal  complexer  is  dan  in  eerste  instan¬ 
ce  aangenomen  werd. 

De  filosofie  van  het  ‘Feasibility-project’  was  indertijd  een  raketsysteem  ongevoe- 
lig  voor  Cook-off  te  maken  door  een  ontsteker  te  bouwen  die  ontsteekt  op  het 
moment  dat  de  hoofdlading  nog  in  een  dusdanige  toestand  verkeerd  dat  een  norma- 
le  verbranding  van  de  motor  waarschijnlijk  is.  Omdat  de  reactiekinetiek  van  de 
stuwstof  niet  voldoende  bekend  was,  was  het  niet  mogelijk  om  gecombineerde 
thermisch-mechanische  berekeningen  uit  te  voeren.  Uit  de  berekeningen  had 
moeten  blijken  of  op  het  moment  van  de  ‘runaway’  van  de  ‘ignitor’,  de  mechani- 
sche  spanningen  in  het  systeem  van  dien  aard  zouden  zijn  dat  een  heftige  reactie 
onwaarschijnlijk  zou  zijn.  Indien  voor  een  vervolg  van  deze  filosofie  gekozen  zou 
zijn,  had  ook  eenvoudig  gebruikgemaakt  kunnen  worden  van  een  ontsteekpil  die 
ontsteekt  bij  lage  temperatuur. 

De  bovenstaande  filosofie  wordt  tegenwoordig  niet  meer  gevolgd  omdat  een 
raketsysteem  dat  bij  130  °C  ontsteekt  niet  gewenst  is  op  een  schip  of  in  een  of 
andere  opslagplaats  voor  munitie-  en  raketsystemen.  Systemen  die  momenteel 
gebruikt  worden  om  de  responsie  van  een  Cook-off  reactie  te  temperen  zijn  onder 
andere  het  gebruik  van  ventgaten  en  andere  technieken  [15, 16, 17,  20, 21, 22]; 
composietmaterialen  [18,  19]  die  bij  voorhoogde  temperatuur  de  opsluiting  ver- 
minderen  en  tenslotte  het  gebruik  van  minder-gevoelige  explosieve  stoffen  en 
stuwstoffen  [23,  24, 25]. 
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Met  het  gebruik  van  de  computermodellen  met  gegevens  uit  de  literatuur,  de 
Tarver-McGuire-reactiekinetiek  en  gecombineerde  drukberekeningen  zijn  goede 
resultaten  bereikt.  Ook  uit  eerdere  berekeningen  is  gebleken  dat  aan  de  hand  van 
dit  model  met  nulde-orde  kinetiek  al  goede  resultaten  behaald  konden  worden  [26]. 
De  temperatuurprofielberekeningen  geven  een  goede  indruk  van  de  temperatuur- 
gradienten  van  het  testitem  die  zeer  goed  overeenkomen  met  de  experimentele 
waarden.  Met  behulp  van  de  drukberekeningen  wordt  een  stopcriterium  voor  de 
simulaties  verkregen  dat  beter  met  de  werkelijkheid  overeenkomt  dan  zuivere 
temperatuurberekeningen . 
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Bijlage  A  ABAQUS  inputfile  en  subroutine  HETVAL 

ABAQUS  INPUTFILE: 


************************************************************************** 
**  TNO-PML  Cook-off  buis 

**  J.H.G.  Scholtes  apr  1997 

**  2 -dim  model  HMX-PBX 

**  aanpassing : integratievolume  via  ABAQUS 

**  aanpassing  cellen  voor  overdracht  tussen  materialen 

**  verschillende  soorten  overdrachten  bij  overgangen 

**  cp  als  functie  van  de  temperatuur  en  Q1=0 . 0418E6 (factor  10 

kleiner  dan  TM) 

************************************************************************** 


♦HEADING 

J.H.G.  SCHOLTES  ' Cook-of f ' buis  (cook5) :Tomg=11.0,Ti=11.0,run  hmx261 
★  * 

♦♦ABAQUS  inputfile  filenaam  co56b2.inp 
** 

***★*★****★*★**★*★★*★★★★ ************************************************** 


**  Temperaturen  in  graden  Celsius 
**  Tijden  in  seconden 
**  Stofgegevens  HMX-PBX  HU  28 -x: 
** 


**  Stofnaam  :  HMXPBX 


**  Tarver  and  McGuire  kinetiek 
♦*  A->B->2C->D 


** 

decomp 

energie 

Q1  : 

-0.0418E6  J/kg 

★  * 

Q2  : 

+1.254E6  J/kg 

** 

Q2  : 

+5.016E6  J/kg 

** 

Aktiv. 

Energie 

Eal: 

26500*R  J/mol 

** 

Ea2 : 

22200*R  J/mol 

** 

Ea3 : 

17200*R  J/mol 

** 

preexp 

factor 

Z1  : 

1 . 412997E21  1/s 

ln(48 . 7 ) 

** 

Z2  : 

1.932160E16  1/s 

ln(37 . 5) 

** 

Z3  : 

1.598361E12  l/s"2 

ln(28 . 1) 

**************** 


**  dichtheid  rho 
**  soort.  warmte  Cp* 
**  soort  geleid  lamb 
**  stralingscoef 
**  epsilon*c 
**  warmteoverdcoef 


1640/1616  kg/m3 

957. 4 +2. 28182 *T 

0.317  J/smK  (T  factor  in  Cp) 

0.0 

4/54E-8  J/sm2K4 

0.5  J/sm2K  (isolation) 


**  explosive  substance  0.87 


***************** 


♦PREPRINT, HISTORY=NO, MODEL=NO 
** 


**  NODES  GENEREREN 

♦NODE 

100,0.0,0.0 

114,0.0170,0.0 

116,0.0175,0.0 

125,0.0280,0.0 

142,0.05,0.0 

900,0.0,0.02 

914,0.0170,0.02 

916,0.0175,0.02 

925,0.0280,0.02 

942,0.05,0.02 
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1000/0.0,0.0205 

1014,0.0170,0.0205 

1016,0.0175,0.0205 

1024,0.0275,0.0205 

1924,0.0275,0.06167 

1925,0.0280,0.06167 

1942,0.0500,0.06167 

7000,0.0,0.27 

7014,0.0170,0.27 

7016,0.0175,0.27 

7024,0.0275,0.27 

*NGEN, NSET=DOPLIOND 

100.114.2 

*NGEN, NSET=DOPREOND 

116.142.2 

*NGEN, NSET=DOPLIMID 

900.914.2 

*NGEN, NSET=DOPREMID 

916.942.2 
*NGEN,NSET=GAPLI 

1000.1014.2 
*NGEN , NSET=GAPRE 

1016.1024.2 

*NGEN, NSET=DOPREBOV 

1924.1942.2 

*NGEN , NSET=BULIBOV 

7000.7014.2 
*NGEN, NSET-BUREBOV 

7016.7024.2 

*NSET , NSET=DOPRRTUS 

925,926,928,930,932,934,936,938,940,942 
*NSET, NSET=DOPRRBOV 

1925,1926,1928,1930,1932,1934,1936,1938,1940,1942 

*NFILL , NSET=BUISES 

GAPLI , BULIBOV , 30,200 

*NFILL, NSET=BUISST 

GAPRE , BUREBOV ,30,200 

*NFILL,NSET=DOPLI 

DOPLIOND, DOPLIMID ,4,200 

*NFILL, NSET=DOPHOEK 

DOPREOND , DOPREMID ,4,200 

♦NFILL, NSET=DOPDRAAD 

DOPRRTOS , DOPRRBOV ,5,200 

*NSET, NSET-NALL 

BUISES , DOPLI , DOPHOEK , DOPDRAAD , BUISST 
*NSET, NSET=HEATER, GENERATE 
3624,7024,200 

*NSET,NSET=NBUISM, GENERATE 

7000.7024.2 

*NSET, NSET=NBUIS5 , GENERATE 

6000.6024.2 

*NSET, NSET=NBUIS4 , GENERATE 

5000.5024.2 

*NSET , NSET=NBOIS3 , GENERATE 

4000.4024.2 

*NSET, NSET=NBUIS2 , GENERATE 

3000.3024.2 

*NSET , NSET=NBUIS1 , GENERATE 

2000.2024.2 

*NSET, NSET=NBUIS0 , GENERATE 

900. 942.2 
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*NSET , NSET=NBUIS00 , GENERATE 

100.142.2 
*NSET, NSET=NPUNT 

NBUISM, NBUIS5 , NBUIS4 , NBUIS3 , NBUIS2 , NBUIS1 , NBUISO , NBUISOO 
*************************** 

**  ELEMENTEN  GENEREREN 

*************************** 

* ELEMENT , TYPE=DCAX4 
101,100,102,302,300 
1001,1000,1002 , 1202,1200 

1017.1016.1018.1218.1216 

926.925.926.1126.1125 

927.926.928.1128.1126 
♦ELGEN, ELSET-EDOND 

101.21.2.2.4.200.200 
♦ELGEN, ELSET=EES 

1001.7.2.2.30.200.200 
*ELGEN, ELSET=EBUIS 

1017.4.2.2.30.200.200 
*ELGEN, ELSET=EDSCHR 

926.1.1.1.5.200.200 
♦ELGEN, ELSET=EDREB 

927.8.2.2.5.200.200 
♦ELEMENT, TYPE=DINTER2A 
923,922,925,1022,1024 
1025,1024,1224,925,1125 
901,900,902,1000,1002 

1015 . 1014 . 1214 . 1016 . 1216 
*ELGEN , ELSET=GAPESS1 

901.11.2.2 

*ELGEN, ELSET=GAPSCHR 

1025.1.2.2.5.200.200 
*ELGEN, ELSET=GAPLENG 

1015.1.2.2.30.200.200 
*ELSET, ELSET=GAPESST 
GAPESS1, 923 
*ELSET,ELSET=GAP 
GAPSCHR , GAPLENG , GAPESST 
*ELSET, ELSET=GAPCOP 
921 

*ELSET, ELSET=GAPSPL 

901,903,905,907,909,911, 913,915,917,919,923 
**★********★**★*★★★*★ 

*ELSET, ELSET=EDOP 
EDOND, EDSCHR, EDREB 
*ELSET, ELSET=ESTAAL 
EDOP, EBUIS 

*ELSET, ELSET=EBUISM, GENERATE 

6801.6823.2 

*ELSET , ELSET=EDRAO , GENERATE 

101.141.2 

★ELSET , ELSET=EDRAR , GENERATE 

141.1741.200 
*ELSET, ELSET=EDBUIT 

1726,1727,1729,1731,1733,1735,1737,1739,1741 
*ELSET , ELSET=HEAT , GENERATE 

2223.6823.200 
*ELSET, ELSET=EALL 
ESTAAL , EES , GAP 

*************************************************************************** 
**  MODEL 

♦PLOT, PLOT  SIZE=1 . 0 , OUTPUT=ASCII , COLORS=4 
COOK-OFF  BUIS  2-D  MODEL  TM-kinetics 
18. ,10. ,13. ,2. 5, 3. 0,2. 0,3. 0,1.0 
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★COLOR  SET , C0L0R=1 , ELSET=ESTAAL 
* COLOR  SET , COLOR=2 , ELSET=EES 
★COLOR  SET , COLOR=4 , ELSET=GAP 
*★  ^SHRINK, FACTOR=0 . 1 

★VIEWPOINT, DEFINITION=MODEL  AXIS  ROTATION 

0.0,0.0,90.0 

★*  *DETAIL, ELSET=CEL 

★DRAW 

************************ 

★★  *DETAIL,ELSET=VLAK 
★*  *DRAW, ELNUM 
**  *SHRINK,FACTOR=0. 0 
**  *VIEWPOINT 
*★  -0.2, -0.5, 1.0,  0.  ,1. ,0. 

**  *DRAW 
** 

********************************************************************** 

★★MATERIAAL  BESCHRIJVING 

★MATERIAL, NAME=HMXPBX 

★CONDUCTIVITY, dependencies=4 

0.31, ,0. ,0. ,0. ,0. 

0.31, ,1000. ,1000. ,1000. ,1000. 

★DENSITY 

1640.0 

******  Solution  dependant  state  variables 

★★★★  1:  Massf raction (A)  is  Y(l)  in  subroutine  HETVAL 

★★★★  2:  Massf raction(B)  is  Y(2)  in  subroutine  HETVAL 

★★★★  3:  Massfraction(C)  is  Y(3)  in  subroutine  HETVAL 

***★  4:  Massfraction(D)  is  Y(4)  in  subroutine  HETVAL 

**★*  5:  Massfraction(A+B+C+D)  is  tot  in  subroutine  HETVAL 
****  6:  pressure  in  integration  area 
★DEPVAR 
6 

★HEAT  GENERATION 
★SPECIFIC  HEAT 
957.4,0.0 
1185.0,100.0 
1414.0,200.0 

*********************************************** 

★MATERIAL , NAME=STAAL 

★CONDUCTIVITY 

50.0 

★DENSITY 

7800 

★SPECIFIC  HEAT 
460.0 
*  * 

** 

★SOLID  SECTION, ELSET=EES,MATERIAL=HMXPBX 

1.0 

** 

★SOLID  SECTION , ELSET=ESTAAL , MATERIAL=STAAL 

1.0 

★INTERFACE, ELSET=GAPLENG 
★GAP  CONDUCTANCE 
3. 65E+02, 0.0005 
★INTERFACE, ELSET=GAPSCHR 
★GAP  CONDUCTANCE 
1000.0 

★INTERFACE, ELSET=GAPCOP 
★GAP  CONDUCTANCE 
1.2E+03 

★INTERFACE, ELSET=GAPSPL 
★GAP  CONDUCTANCE 
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*** 

************************************************************************* 
**  KINETISCHE  VERGELI JKINGEN 

★USER  SUBROUTINE 

INCLUDE  1 /local/usrl/scholtes/druk/tmco56b2 . f ' 
********************************************************* 

**  INITIELE  CONDITIES 

★INITIAL  CONDITIONS , TYPE=FIELD , VARIABLES 

NALL, 0 . 

★INITIAL  CONDITIONS, TYPE=TEMPERATURE 
NALL, 11.0 

★INITIAL  CONDITIONS, TYPE=SOLUTION 
EES, 1.0, 0.0, 0.0, 0 .0, 1. 0 

★AMPLITUDE, NAME=LINEAR,DEFINITION=TABULAR, TIME=TOTAL  TIME, VALUE= ABSOLUTE 
0.0,11.0,25000.0,255.75 

***********************  =  T=ll . 0+0 . 00979 *t 
★RESTART , WRITE , FREQUENCY'S 
** 

**************************************************  STEP  gedeelte 
★★  begin  voor  integratie  volume  ******** 

★*  nieuw  ivm.  het  integratie  volume 

★★  een  zeroe  step  doet  niets,  Dit  is  nodig  omdat  IVOL  via  getvrm 
**  na  increment  #1  doorgegeven  wordt 
*  * 

★step 

★heat  transfer,  steady  state 
★boundary, AMPLITUDE=LINEAR 
HEATER, 11,  ,0.0 
★end  step 

*******************  einde  voor  int  volume*************** 

**  INC=number  of  increments; 

** 

★STEP, INC=250 

**  DELTMX=maximum  temperature  change  in  time- step 
**  Next  line:  Initial  step-size  and  Time-limit: 

★HEAT  TRANSFER , DELTMX=4 . 0 , END=PERIOD 
1.0,25000.0, , , 

**  time-limit  is  25000  SECONDE  ** 

************************************************* 

★BOUNDARY , AMPL I TUDE=L I NEAR 
HEATER, 11,  ,0.0 
*  * 

**  Ambient  temperature  and  Heat  Transfer  Coefficient  en  straling  : 

★FILM 

EDRAO , Fl , 11 .0,0.5 
EDRAR, F2, 11 . 0, 0 . 5 
EDBUIT,F3, 11.0, 0.5 
★RADIATE, ZERO=- 27 3 .16 
EDRAO, Rl, 11.0,4. 54E- 18 
EDRAR, R2, 11.0,4 . 54E-18 
EDBUIT, R3 ,11.0,4. 54E- 18 

************************************************ 

**  invoeren  loads  om  bij  mater iaal  berek  ufield  aan  te  roepen 
************************************************ 

★field, user, variable=l 
NALL 

★field, user , variable=2 
NALL 

★field, user , variable=3 
NALL 

★field, user, variable=4 
NALL 

************************************************* 
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★PRINT, FREQUENCY=1 

★NODE  PRINT , NSET-NALL , FREQUENCY=3 

NT 

★NODE  FILE, NSET=NALL, FREQUENCY=1 
NT 

**  *EL  PRINT 
**  I VOL 

★★★PLOT, PLOT  SIZE=1 . 0 , FREQUENCY=20 , OUTPUT=ASCII 
★★TEMPERATUUR  CONTOURPLOT 
**30. 0,10. 0,26. 0,5. 0,3. 0,3. 0,3. 0,1.0 
★* *VIEWPOINT , DEFINITION=MODEL  AXIS  ROTATION 
*★0.0, 0.0, 90.0 

★★★PLOT  MODE , ELEMENT=NO , NODE=NO , FILL=NO 

★★★CONTOUR 

★★TEMP 

**  *EL  FILE, FREQ=1, ELSET= 

**  TEMP 
★ENDSTEP 

************************************************** 
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SUBROUTINE  HETVAL: 

c* ******************************************************************* 
C  subroutine  tmco56b2.f  for  TARVER -McGuire  kinetics  for  HMXPBX  * 

c  subroutinte  voor  drukberekening  voor  co56b2.inp  met  correctie  * 

c  in  de  REP-formule:  ->  "+"  en  versie  5.6  van  ABAQUS  * 

C  Verdere  aanpassing  voor  Gap  conductance  uitvoering  * 

C******************************************************************** 

SUBROUTINE  HETVAL ( CMNAME , TEMP , TIME , DTIME / SVAR, FLUX, PREDEF ) 

C  TEMP ( 1 )  is  temperature  after  this  increment 

C  TEMP (2)  is  temperature  increment 

C  TIME ( 1 )  is  the  step  time  at  the  end  of  this  increment 

C  TIME ( 2 )  is  the  total  time  after  this  increment 

C  DTIME  is  the  time  step  dt 
C  kinetics  by  Tarver  en  McGuire  model 
IMPLICIT  REAL *8(A-H,0-Z) 

CHARACTER *8  CMNAME 

DIMENSION  TEMP ( 2 ) , TIME ( 2 ) , SVAR ( 6 ) , FLUX ( 2 ) , PREDEF ( * ) 

INTEGER  kmax,  kount ,  nbad,  nok,  ILEL,  ILNO,  IDONE,  kinds 
REAL* 8  dxsav,eps,hstart,xl,x2/y(4) , xp, yp, T, tot , 

*el , e2 , e3 , zl , z  2 , z3 , rql , rq2 , rq3 , RP , FTOT , RUITK (3,350) 

& , RPVT1 , RPVT2 , RVTl , RVT2 , RPGLl , RPGL2 , RVOL 
INTEGER  KMAXX, NMAX 
PARAMETER  (KMAXX=200 , NMAX=50 ) 

PARAMETER (el=26 5 00 . ,e2=22200. , e3=17200 . , zl=l . 412997D21, 

*z2=l . 93216D16 , z3=l . 598361D12) 

COMMON  /KPATH/  kmax, kount, dxs a v, xp (KMAXX) , yp (NMAX, KMAXX) 

COMMON  /KTEMP/  T 

common/celgi/iduml (7 ) , iedbr, idum2 ( 134 ) 
common/cmatgi/idum3 ( 23 ) , imtemp, idum4(66) 
common/count/kinc , idum5 (58) 
c 

common/cedat/idum6 ( 2 ) , ienod 
common/cns/idum7 ( 1 ) , incrd , inu , idum8 ( 1 ) 
common/cuvarm/keltemp, kintuv, ksptuv 

COMMON  /KGPRES/RPABl , RPAB2 , RPVTl , RPVT2 , RVTl , RVT2 , ILEL , ILNO 
COMMON  /KUITK/RUITK,  IDONE,  kinds 
c****************************************************** 
c  ****  aanpassing  voor  integratie  volume  via  abaqus  *** 
c  common/cmats/idumx ( 4 ) , xdum , j  dum ( 59 ) , luvarm 

c  ***  aanpassing  voor  abaqus  5.5.  (59)  moet  (60)  zijn! ! 

C****************************************************** 

common/cmats/idumx ( 4 ) , xdum , j  dum (60), luvarm 
character*8  flgray(15) 
dimension  arraygs (15) , jarray(15) 

C  , idumx(4 ) , j dum (60 ) 

c  **********einde  aanpassing  IVOL 
c 

c  ************************************************************** 

c  jelno  :  current  element  number 
c  nnod  -  no.  of  nodes  on  the  element 
c  ocoord(2)  =  original  nodal  coords 
c  nnum  =  current  node  no.  (user  defined) 

c  inno  =  node  numbers  of  the  element 

c  kinc  =  current  increment  number 

0  ************************************************************** 

c  character *8  cmname 

dimension  coords (3) 

dimension  array (10000) ,ocoord(2) , fcoord(2) ,disp(2) 

C  dimension  temp(2) ,svar(l) ,predef (*) ,time(2) , flux(l) 

dimension  nodes (27) 

C 

EXTERNAL  stifbs , derivs 

C**************************************************************** 
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C  startwaarden  van  toestandsvariabelen  en  startwaarden 

C  voor  het  stelsel  diffvlgn  toekennnen 

C  en 

C  initial  massf ractions  in:  *INITIAL  CONDITIONS, TYPE=SOLUTION 

C**************************************************************** 

eps=lD-3 

hstart=DTIME 

kmax=0 

y ( 1 ) =SVAR ( 1 ) 
y ( 2 ) =SVAR ( 2 ) 
y (3 )=SVAR( 3 ) 
y(4 )=SVAR(4 ) 

T=TEMP ( 1 ) - 0 . 5  *  TEMP (2)+273.15 

xl=TIME(l)-DTIME 

x2=TIME (1) 

C  write ( 8 , * )  T, y (1) , y(2 ) , tot, xl, x2 

C************************************************************** 

C  aanroep  van  het  oplossen  van  diff  vgl  van  reactiekinetiek 

C************************************************************** 

call  odeint ( y , 4 , xl , x2 , eps , hstart , 0 . , nok , nbad , derivs , stif bs ) 
tot=y ( 1 ) +y ( 2 ) +y ( 3 ) +y ( 4 ) 

SVAR  ( 1 )  ==y  ( 1 ) 

SVAR ( 2 ) =y ( 2 ) 

SVAR ( 3 ) =y ( 3 ) 

SVAR ( 4 ) =y ( 4 ) 

SVAR ( 5 ) =tot 

C*****  berekenen  van  de  vrijgekomen  warmte  in  eenheidscel 
rql=-0 . 0418E6*zl*dexp( -el/T) *y (1) 
rq2=+l . 254E6*z2*dexp( -e2/T) *y ( 2 ) 
rq3=+5. 016E6*z3*dexp(-e3/T)*(y(3)**2) 

FLUX (1)=0. 87*1640* (rql+rq2+rq3 ) 

FLUX ( 2 ) =0 . 87*1640* (1/(T*T) ) * (el*rql+e2*rq2+e3*rq3 ) 

C  write(8 , * )  FLUX(l) ,FLUX(2) ,y(l) ,y(2) ,y(3) ,y(4) , tot,x2 

C************************************* 

C*  berekeningen  voor  node  en  elementen 
C* ************************************ 

c 

ipoint 1= - ( imtemp - iedbr ) 
c 

call  acopdj (predef (ipointl+1) , jelno) 
c 

ipoint2=- ( imtemp- iedbr )+ienod-l 
c  nnod=predef (ipointl+5) 

call  acopdj (predef ( ipointl+5 ) , nnod) 
do  i=l,nnod 

call  acopdj (predef ( ipoint2+i-l) , inno) 
call  dbnd(0, array , inno, ircd) 
call  acopdi(array(2) , nnum, 1) 
nodes ( i ) =nnum 
ocoord ( 1 ) =array ( incrd ) 
ocoord ( 2 ) =array ( incrd+1 ) 
c  include  the  following  line  for  a  3-d  model 
c  ocoord ( 3 )=array ( incrd+2) 

c  write(7, *)kinc, jelno, nnum, ocoord(l) ,ocoord(2) 

enddo 

C+++ 

C  Here  is  what  you’ll  find  in  PREDEF 
C 

C  PREDEF (1)  =  Current  node  #.  Since  PREDEF  is  a  REAL* 8  array  you 

C  should  convert  it  to  an  integer.  I  used  the  built-in 

C  function  DINT  for  doing  it  and  written  the  result  to 

C  nodnum . 

C 

C  PREDEF ( 2 )  =  X- coordinate  of  this  node 
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C  PREDEF ( 3 )  =  Y- coordinate  " 

C  PREDEF ( 4 )  =  Z- coordinate  " 

C 

c--- 

nodnum=dint ( predef ( 1 ) ) 

C* **************************************** * 

C  ***  verdere  aanpassingen  voor  IVOL******* 

C  ***  ophalen  van  volume  van  eenheidscel 
C* ***************************************** 


luvarm=l 

call  getvrm( ' IVOL' , arraygs , j array, flgray, jrcd) 

luvarm=0 

uvar=arraygs ( 1 ) 

C  ********einde  IVOL  aanpassing  ********** 

C 

C  Here  some  write's  to  check  that  the  sub  does  the  right  things... 

C  You'll  find  the  output  in  the  .msg  file  (Unit  7) 

C - 


c 

write (7 , * ) 

* - From 

HETVAL . 

c 

write ( 7 , * ) 

V  f 

c 

write(7, *) 

' Increment# : 

' , kinc 

c 

write (7 , * ) 

'Accumulated  time: 

' , time(2) 

c 

write (7 , * ) 

'El#: 

' , jelno 

c 

write(7, *) 

' Nodes  forming  El . : 

' , ( nodes ( i ) , i=l , nnod ) 

c 

write(7, *) 

'  ' 

c 

write (7, *) 

»  i 

c 

write (7, *) 

'Volume  at  iNt.pnt. : 

' , uvar 

c 

write(7, *) 

'Actual  node#: 

' , nodnum 

c 

write ( 7 , * ) 

'Coordinates: 

' , predef ( 2 ) , predef ( 3 ) 

c 

write (7 , *) 

'  ' 

c 

write (7 , *) 

'  ++  +  +  +++++  +++++++++-4 

-++++++++++++++++++++ ' 

c  calculation  of  the  pressure  in  node-area 

C  RELEASED  GASSES  IN  THE  2ND  AND  3TH  REACTIONSTEP 


C****  met  behulp  van  de  gasproductief actor  de  druk  berekenen 
C****  in  deze  eenheidscel 

q* ********************************************************* 


FTOT=SVAR ( 3 ) +SVAR ( 4 ) 

CALL  KPRESS ( FTOT, 0 . , TEMP ( 1 ) , RPGl , RPAB1 , RPG2 , RPAB2 ) 
c  WRITE  (7,*)  "PRESSURE",  RPGl 

C* ********************************************************* 


C  *  *  *  * * VOORDAT  ITERATIE  OPNIEUW  BEGINT  HET  PRINTEN  VAN  DE  GLOBALE  DRUK  *** 


C *********************************************************  * 

C  IF  ( (jelno.EQ. 6813) .and. (nodnum.EQ. 7014) )  then 

C . IF  (kinc  .GT  .kinds)  then 

RPGL1=RPVT1/RVT1 

RPGL2=RPVT2/RVT2 

WRITE ( 7  ,  *  )  '  ***  Pgl=  ' ,RPGLl,RPGL2 

RUITK ( 1 , kinc ) =TIME ( 2 ) 

RUITK ( 2 , kinc ) =RPGL1 

RUITK ( 3 , kinc ) =RPGL2 

IF  (IDONE.EQ.l)  GOTO  99 

IF  (RPGL1.GT. 400.0)  THEN 

WRITE ( 7 , * )  ' Pglobaal  druk-tijd: ' 

write (7 , 101)  (RUITK ( 1 , IT) , RUITK ( 2 , IT) , RUITK ( 3 , IT) , IT=l,kinc) 
WRITE ( 7 , * )  'EINDE  UITVOERTABEL ' 

IDONE=l 

ENDIF 

99  ILEL=jelno 
ILNO=nodnum 
RPVT1=0 . 

RPVT2=0. 

RVT1=0 . 

RVT2=0 . 
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END  IF 

Q*  ******************************************************  * 

C****  ophogen  van  de  som  voor  de  globale  drukberekening 
q*  ******************************************************** 

RVOL=uvar 

RPVTl=RPVTl+RPGl*RVOL*RPABl/( TEMP ( 1 ) +273 . 15 ) 
RPVT2=RPVT2+RPGl*RVOL*RPAB2/ ( TEMP( 1 )+27  3.15) 

RVT1=RVT1+RV0L*RPAB1/ ( TEMP ( 1 ) +27  3 . 15 ) 

RVT2  =RVT  2  +RVOL *  RP AB  2/ ( TEMP ( 1 ) +  2  7  3 . 1 5 ) 
kincls=kinc 

Q* **********************************************  ********* 

100  FORMAT (X/ 216 , 2F9 . 6 , 4F10 .2) 

101  FORMAT (X , 3F9 . 3) 

C  **  indien  de  cel  tot  de  centrumrij  behoort  uitprinten: 

C* ******************************************************* 

IF( jelno.LT. 6800)  goto  200 

WRITE  (7,100)  j elno , kinc , SVAR ( 3 ) , SVAR ( 4  )  , TEMP ( 1 ) , TIME ( 2 ) , 

&RPG1 , RPG2 
200  RETURN 
END 

Q* *********************************************  ****** 

c  Niet  aankomen  via  abaqus  *** 
c************************************************ 

C  UFIELD 

c - 

subroutine  ufield(f ield, kfield, nsecpt, kstep, kinc, time, node, 

& coords , temp , dtemp ) 
c  include  ' aba_param . inc ' 

IMPLICIT  REAL*8(A-H,0-Z) 

dimension  field (nsecpt) , time(2) , coords (3) , temp(nsecpt) 
dimension  dtemp (nsecpt) 

C+++ 

C 

C  if  we  are  processing  variable=l  from  the  *FIELD  then  write  the  node 

C  number  to  field(l) .  This  will  be  written  subsequently  to  PREDEF(l) 

C 

C-  - 

if  (kfield  .eq.  1)  then 
f ield(l)=node*l . 
else 

C+++ 

C  if  we  are  processing  variable=2,3  or  4  then  write  X,  Y  or  Z  to 

C  field(l).  We'll  find  these  values  in  PREDEF ( 2 ) ,  PREDEF ( 3 )  and 

PREDEF ( 4 ) 

C  within  HETVAL 

C  — 

field ( 1 ) =coords ( kf ield- 1 ) 
endif 

return 

end 

C***** *****************  ********  ********************** 

SUBROUTINE  KPRESS ( RF , REPSIL , RT , RPG1 , RPABl , RPG2 , RPAB2 ) 

C  SUBROUTINE  PRESSURE  ESTIMATION 

C  DEZE  SUBROUTINE  MOET  DE  DRUK  IN  EEN  CEL  BEREKENEN 

C  VIA  DE  METHODE  VAN  BKW.  DE  METHODE  WORDT  GEBRUIKT  VOOR 

C  HET  COMPUTERPROGRAMMA  TIGER  EN  DE  INVOERWAARDEN  KOMEN 

C  HOOFDZAKELIJK  UIT  HET  ARTIKEL: 

C  "CONSTITUTIVE  MECHANICAL  MODEL  FOR  ENERGETIC  MATERIALS" 
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C  VAN  M.L  HOBBS,  MELVIN  R.  BAER  AND  R.J.  GROSS 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL *8  RALFA, RBETA, R KAPPA, RTHETA, RCV, RXC, RHOCOl, RHOTH, RMW, RYKI , 
&RR , RHOGO , RHOC , RHOG , RPABO , RF , RLAMBA , RPSI , RN , RKT , RY , RPGO , RPG1 , RPG2 
&RX, RZ , RBO , RB , RAO , RA, REPSIL , RT , RDC , RFE , RATBB , RPAB1 , RPAB2 , RHOC 02 
PARAMETER  (RALFA=0 . 5 , RBETA=0 . 284 , RKAPPA=11 . 85 , RTHETA=10200 . , 
&RDC=1 . 1, RCV=1 . 0032, RXC=0 . 32, RHOTH=l . 64 , RMW=33 . 3 , 

&RYKI=471 . , RR=8 . 3144 , RHOGO=1 . 29D- 03 , RN=8 . 3 , RKT=14100 . , RY=100 . , 
&RPG0=0 . 1 , RHOC01=1 . 6166 , RHOC02=1 . 6266 ) 

C  INITIAL ISATIE  VAN  PARAMETERS  met  RHOCOl 

C  EFFECT I EVE  OMZETGRAAD : 

RFE=RF*(1. -RXC) 

RPAB0=1 . “ ( RHOCO 1/RHOTH ) 

RLAMBDA= (RHOCOl/RHOGO ) * ( (1 ./RPABO) -1 . )+1.0 
RPSI=1 . +RFE* (RLAMBDA-1 . ) 

C  RATIO  =RATBB=B0**3/B**3 . 

RATBB-1 . / ( 1 . +REPSIL ) *  *3 . 

RPABl=RATBB* (RPABO+RFE* ( 1 . -RPABO ) ) 

C  WRITE (7,*)  RFE, RPABO , RLAMBDA, RPSI, RATBB, RPAB 

C 

C  DE  HOOFDVERGELI JKINGEN  OM  DE  DRUK  TE  BEREKENEN : 

C 

RHOG=RHOGO*RPABO*RPSI/(RPAB1* ( 1 . +REPSIL) ) 
RX=RHOG*RKAPPA*RYKI/(RMW*(273 . 15+RT+RTHETA) **RALFA) 

RZ  =  1 . +RX*DEXP ( RBETA*RX ) 

RPGl=RZ*RHOG*RR* (RT+273 . 15)/RMW 

RHOC=RHOG* ( ( RLAMBDA/RPS I ) - 1 . )/( (1 ./RPABl) -1 . ) 

RPC=RHOC* RDC* RCV* (RT+273 . 15 ) +RKT* ( (RHOC/RHOCOl ) **RN- 1 . 0 )/RN 
C  WRITE ( 7 , * )  RHOG , RX , RZ , RPG, RHOC , RPC 

q* * ***************************************** * 

C  berekeningen  van  parameters  van  de  hoofdvgln 

Q  +  *** *****************************  ************ 

C  EFFECTIEVE  OMZETGRAAD : (  MET  RHOC02) 

RFE=RF*(1. -RXC) 

RPAB0=1 . - ( RHOCO 2/RHOTH) 

RLAMBDA= ( RHOCO 2/RHOG0)*( (1. /RPABO) -1. )+1.0 
RPSI=1.+RFE*( RLAMBDA-1. ) 

C  RATIO  =RATBB=B0  *  *  3/B  *  *  3 . 

RATBB=1 . /( 1 . +REPSIL) *  *3 . 

RPAB  2=RATBB  * (RPABO+RFE* ( 1 . -RPABO ) ) 

C  WRITE (7,*)  RFE, RPABO, RLAMBDA, RPSI, RATBB, RPAB 

Q** *****************************************  ***************** 

C  DE  HOOFDVERGELI JKINGEN  OM  DE  DRUK  TE  BEREKENEN:******** 

q*** ***********************************  ********************** 
RHOG=RHOGO  *RP AB  0  *RPS 1/ ( RPAB  2  * ( 1 . +REPSIL ) ) 

RX=RHOG*RKAPPA*RYKI/ ( RMW* (273. 15+RT+RTHETA ) *  * RALFA ) 

RZ=1 . +RX*DEXP ( RBETA*RX ) 

RPG2=RZ*RHOG*RR* (RT+273 . 15 )/RMW 

RHOC=RHOG*( (RLAMBDA/RPSI)-!. )/( (l./RPAB2)-l. ) 

RPC=RHOC*RDC*RCV* ( RT+273 . 15 ) +RKT* ( (RHOC/RHOC02 ) *  *RN- 1 . 0 ) /RN 
C 

RETURN 

END 

Q*** ****************************************  ********* 

q***  oude  subroutine  om  volume  van  cel  te  bepalen 
q* ************************************  *************** 

SUBROUTINE  KVOLUM ( JELNO , NODNUM , RVOL ) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  JELNO, NODNUM, IH, IK 
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REAL *8  RVOL,RAR(7,4) 

C 

RAR( 1,1)=1.6597D"08 

RAR(1, 2 ) =6 .1942D-08 

RAR(1, 3 )=1 . 6597D-08 

RAR(l,4)-6.1942D-08 

RAR (2,l)=9.5137D-08 

RAR(2, 2)=1.4048D-07 

RAR(2,3)=9.5137D-08 

RAR (2/ 4 )=1 . 4048D-07 

RAR (3,l)=1.7368D-07 

RAR(3,2)=2.1902D-07 

RAR (3,3) =1 . 7  368D- 07 

RAR (3,4)=2.1902D-07 

RAR (4 ,l)=2.5222D-07 

RAR(4,2)=2.9756D-07 

RAR (4 / 3 )—2 . 5222D-07 

RAR(4, 4) =2. 975 6D -07 

RAR(5, l)=3.3076D-07 

RAR (5,2)=3.7610D-07 

RAR( 5 , 3 ) =3 . 3076D- 07 

RAR(5,4)=3 .7610D-07 

RAR(6, l)=4.0930D-07 

RAR ( 6 , 2 ) =4  . 5464D- 07 

RAR ( 6 / 3 )=4 . 0930D- 07 

RAR(6,4)=4 . 5464D-07 

RAR (7,1) =4 . 8784D-07 

RAR (7 / 2 )=5 . 3318D-07 

RAR (7/3)=4.8784D~07 

RAR(7, 4)=5 .3318D-07 

IH-JELNO- ( JELNO/100 ) *100 

IF  ( (JELNO-NODNUM) .EQ.l)  THEN 

IK=1 

END  IF 

IF  ( (JELNO-NODNUM) . EQ . 0 )  THEN 

IK=2 

END  IF 

IF  ( (JELNO-NODNUM) .EQ. 99)  THEN 

IK-3 

END  IF 

IF  ( (JELNO-NODNUM) .EQ. 100)  THEN 

IK-4 

ENDIF 

RVOL—RAR ( IH / IK) 

RETURN 

END 

C* *************************************************** 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  sl#  +  ] . .  . 

C  Used  subroutines 

SUBROUTINE  odeint ( y start , nvar , xl , x2 , eps , hi , hmin , nok , nbad , der ivs , 
*rkqs) 

IMPLICIT  REAL* 8  (A-H,0"Z) 

INTEGER  nbad , nok , nvar , KMAXX , MAXSTP , NMAX 
REAL* 8  eps , hi / hmin, xl , x2 , ystart (nvar ) , TINY , T 
COMMON  /KTEMP/  T 
EXTERNAL  derivs , rkqs 

PARAMETER  (MAXSTP-10000 , NMAX-50 , KMAXX-200 , TINY-1 . D- 30 ) 

INTEGER  i, kmax, kount, nstp 

REAL* 8  dxsav, h,hdid, hnext , x, xsav, dydx(NMAX) ,xp(KMAXX) , y ( NMAX ) , 
*yp ( NMAX , KMAXX )  , yscal ( NMAX ) 

COMMON  /KPATH/  kmax , kount , dxsav, xp, yp 
x=xl 

h=dsign(hl, x2-xl) 
nok=0 
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nbad=0 

kount=0 

do  11  i=l,nvar 
y(i)=ystart(i) 

11  continue 

if  (kmax.gt.O)  xsav=x-2 . *dxsav 
do  16  nstp=l / MAXSTP 

call  derivs (x, y, T, dydx) 
do  12  i=l,nvar 

yscal (i)=dabs(y(i) )+dabs(h*dydx(i) )+TINY 

12  continue 

if ( kmax . gt . 0 ) then 

if (dabs (x-xsav) .gt .dabs (dxsav) )  then 
if ( kount . It . kmax- 1 ) then 
kount=kount+l 
xp( kount )=x 
do  13  i=l,nvar 
yp(i, kount )=y(i) 

13  continue 
xsav=x 

endif 

endif 

endif 

if ( (x+h-x2 ) * (x+h-xl) . gt . 0 . )  h=x2-x 

call  rkqs ( y , dydx , nvar , x, h , eps , yscal , hdid , hnext , derivs ) 
if ( hdid . eq . h ) then 
nok=nok+l 
else 

nbad=nbad+l 

endif 

if ( (x-x2) *(x2-xl) . ge. 0 . )then 
do  14  i=l,nvar 
ystart(i)=y(i) 

14  continue 

if ( kmax . ne . 0 ) then 
kount=kount+l 
xp( kount )=x 
do  15  i=l,nvar 
yp  ( i , kount ) =y ( i ) 

15  continue 
endif 
return 

endif 

if (dabs (hnext) . It .hmin)  pause 
*'stepsize  smaller  than  minimum  in  odeint' 
h=hnext 

16  continue 

pause  ’too  many  steps  in  odeint' 

return 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  sl#+] . . 

C  ******************* 

SUBROUTINE  st if bs (y , dydx , nv, x, htry , eps , yscal , hdid, hnext, derivs ) 
IMPLICIT  REAL *8  (A-H,0-Z) 

INTEGER  nv , NMAX , KMAXX , IMAX 

REAL* 8  eps , hdid , hnext , ht ry , x , dydx ( nv ), y ( nv ), yscal ( nv ) , SAFEl , SAFE2 , 
*REDMAX , REDMIN , TINY , SCALMX , T 
COMMON  /KTEMP/T 
EXTERNAL  derivs 

PARAMETER  (NMAX=5 0 , KMAXX-7 , IMAX=KMAXX+1 , SAFEl= . 25 , SAFE2= . 7 , 
*REDMAX=1 . e- 5 , REDMIN= . 7 , TINY=1 . e- 30 , SCALMX= . 1 ) 
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CU  USES  derivs , jacobn, simpr, pzextr 

INTEGER  i, iq, k, kk, km, kmax, kopt, nvold, nseq( IMAX) 

REAL* 8  epsl, epsold,  errmax,  fact,h,  red, scale, work , wrkmin,xest,xnew, 
*a ( IMAX) , alf (KMAXX, KMAXX) , dfdx(NMAX) , dfdy (NMAX, NMAX) , err (KMAXX) , 
*yerr (NMAX) , ysav(NMAX) , yseq (NMAX) 

LOGICAL  first, reduct 

SAVE  a , alf , epsold , f irst , kmax , kopt , nseq , nvold , xnew 
DATA  first/. true ./, epsold/- 1 ./, nvold/- 1/ 

DATA  nseq  /2 , 6 , 10 , 14 , 22 , 34 , 50 , 70/ 
if ( eps . ne . epsold . or . nv . ne . nvold ) then 
hnext=-l . e29 
xnew=-l.e29 
epsl=SAFEl*eps 
a(l)=nseq(l)+l 
do  11  k=l , KMAXX 

a(k+l)=a(k)+nseq(k+l) 

11  continue 

do  13  iq=2, KMAXX 
do  12  k=l, iq-1 

alf (k, iq)=epsl**( (a(k+l) -a(iq+l) )/( (a(iq+l) -a(l)+l. )*(2*k+ 

*1))) 

12  continue 

13  continue 
epsold=eps 
nvold=nv 
a(l)=nv+a (1) 

do  14  k=l, KMAXX 

a(k+l)=a(k)+nseq(k+l) 

14  continue 

do  15  kopt=2 , KMAXX -1 

if ( a ( kopt+1 ) . gt . a ( kopt ) *alf ( kopt - 1 , kopt ) ) goto  1 

15  continue 

1  kmax=kopt 
endif 
h=htry 

do  16  i=l,nv 
ysav(i)=y(i) 

16  continue 

call  jacobn (x, y, T, dfdx, dfdy, nv, nmax) 
if (h . ne . hnext . or . x . ne . xnew) then 
f irst= .true. 
kopt=kmax 
endif 

reduct= . false . 

2  do  18  k=l,kmax 

xnew=x+h 

if (xnew. eq .x) pause  'stepsize  underflow  in  stifbs’ 
call  simpr ( ysav , T , dydx , dfdx , dfdy , nmax , nv, x , h , 

*nseq ( k ) , yseq , derivs ) 
xest= ( h/nseq ( k ) ) *  *  2 
call  pzextr ( k , xest , yseq , y , yerr , nv, T ) 
if (k.ne.l)then 
errmax=TINY 
do  17  i=l,nv 

errmax=max ( errmax , abs ( yerr ( i ) /yscal ( i ) ) ) 

17  continue 
errmax=errmax/eps 
km=k-l 

err ( km )  =  ( errmax/SAFEl ) *  * ( 1 . / ( 2  *km+l ) ) 
endif 

if ( k . ne . 1 . and . ( k . ge . kopt - 1 . or . f irst ) ) then 
if (errmax . It . 1 . )goto  4 
if ( k . eq . kmax . or . k . eq . kopt+1 ) then 
red=SAFE2/err ( km ) 
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goto  3 

else  If ( k . eq . kopt ) then 

if (alf (kopt-1, kopt) . It .err (km) )then 
red=l ./err (km) 
goto  3 
endif 

else  if (kopt. eq.kmax) then 

if (alf (km, kmax-1) . It . err (km) ) then 
red=alf (km, kmax-1) *SAFE2/err (km) 
goto  3 
endif 

else  if (alf (km, kopt) .lt.err(km) )then 
red=alf (km, kopt-l)/err (km) 
goto  3 
endif 
endif 

18  continue 

3  red=min(red, REDMIN) 
red=max ( red , REDMAX ) 
h=h*red 

reduct = . true, 
goto  2 

4  x=xnew 
hdid=h 

f irst= . false . 
wrkmin=l . e35 
do  19  kk=l,km 

fact=max(err (kk) , SCALMX) 
work=f act*a ( kk+1 ) 
if (work . It . wrkmin ) then 
scale=fact 
wrkmin=work 
kopt=kk+l 
endif 

19  continue 
hnext-h/scale 

if ( kopt . ge . k . and . kopt . ne . kmax . and . . not . reduct ) then 
fact=dmaxl(scale/alf (kopt-1, kopt) , SCALMX) 
if (a (kopt+1) *fact . le .wrkmin) then 
hnext=h/fact 
kopt=kopt+l 
endif 
endif 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0-,. 

SUBROUTINE  simpr (y, T, dydx, dfdx, df dy , nmax, n,xs,htot, nstep,yout, 
*derivs ) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  n, nmax, nstep, NMAXX 

REAL* 8  htot,xs,dfdx(n) , dfdy( nmax, nmax) ,dydx(n) ,y(n) ,yout(n) ,T 
EXTERNAL  derivs 
PARAMETER  (NMAXX=50) 

CU  USES  derivs , lubksb, ludcmp 
INTEGER  i, j ,nn, indx( NMAXX) 

REAL* 8  d,h,X,a(NMAXX, NMAXX) , del (NMAXX ), ytemp ( NMAXX ) 
h=htot/nstep 
do  12  i=l,n 
do  11  j=l,n 

a(i, j )=-h*dfdy(i, j ) 

11  continue 

a( i, i)=a (i, i)+l . 

12  continue 

call  ludcmp ( a, n, NMAXX, indx, d) 
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do  13  1=1, n 

yout(i)=h*(dydx(i)+h*dfdx(i) ) 

13  continue 

call  lubksb ( a , n , NMAXX , indx , yout ) 
do  14  i=l,n 

del(i)=yout(i) 
ytemp ( i ) =y ( i ) +del ( i ) 

14  continue 
x=xs+h 

call  deri vs (x, ytemp, T, yout ) 
do  17  nn=2,nstep 
do  15  i=l,n 

yout ( i ) =h*yout ( i ) -del ( i ) 

15  continue 

call  lubksb ( a, n, NMAXX, indx, yout) 
do  16  i=l, n 

del ( i ) =del ( i ) +2 . *yout ( i ) 
ytemp ( i ) =ytemp ( i ) +del ( i ) 

16  continue 
x=x+h 

call  der i vs (x, ytemp, T, yout) 

17  continue 

do  18  i*l,n 

yout ( i ) =h*yout ( i ) -del ( i ) 

18  continue 

call  lubksb ( a, n, NMAXX, indx, yout) 
do  19  i=l,n 

yout ( i ) =y temp ( i ) +yout ( i ) 

19  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0-,. 
SUBROUTINE  pzextr ( iest , xest , yest , yz , dy , nv, T ) 

IMPLICIT  REAL  *  8  (A-H,0-Z) 

INTEGER  iest,nv, IMAX,NMAX 

REAL* 8  xest, dy(nv) ,yest(nv) ,yz(nv) ,T 

PARAMETER  ( IMAX=13 , NMAX=50 ) 

INTEGER  j , kl 

REAL *8  delta , f 1 , f 2 , q , d ( NMAX ) , qcol ( NMAX , IMAX ) , x ( IMAX ) 
SAVE  qcol,x 
x(iest)=xest 
do  11  j=l,nv 
dy(j)=yest(j) 
yz( j )=yest( j ) 

11  continue 

if (iest.eq.l)  then 
do  12  j=l,nv 

qcol ( j , 1 ) =yest ( j ) 

12  continue 
else 

do  13  j=l,nv 
d(j)=yest( j) 

13  continue 

do  15  kl=l, iest-1 

delta=l ./(x( iest-kl) -xest) 
f l=xest*delta 
f 2=x( iest-kl) *delta 
do  14  j=l,nv 
q=qcol( j , kl) 
qcol ( j , kl ) =dy ( j ) 
delta=d( j ) -q 
dy ( j )=f l*delta 
d( j )=f 2*delta 
yz( j )=yz( j )+dy( j ) 
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14 

continue 

15 

continue 

do  16  j=l,nv 

qcol( j, iest)=dy( j ) 

16 

continue 

endif 

return 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  0-,. 

C  ******************* 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  sl#+] . . . 
SUBROUTINE  lubksb (a, n, np, indx, b) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  n,np,indx(n) 

REAL* 8  a(np,np) ,b(n) 

INTEGER  i, ii, j , 11 
REAL *8  sum 
ii-0 

do  12  i=l,n 
ll=indx( i) 
sum=b(ll) 
b(ll)=b(i) 
if  (ii.ne. 0)then 
do  11  j=ii,i-l 

sum=sum-a(i, j )*b( j ) 

11  continue 

else  if  (sum.ne.O.)  then 
ii=i 
endif 
b(i)=sum 

12  continue 

do  14  i=n, 1, -1 
sum=b(i) 
do  13  j=i+l,n 

sum=sum-a(i, j )*b( j ) 

13  continue 

b( i)=sum/a( i, i) 

14  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  sl#+] . . . 
SUBROUTINE  ludcmp (a, n , np, indx, d) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  n , np , indx ( n ) , NMAX 

REAL *8  d ,  a  ( np , np ) , TINY 

PARAMETER  (NMAX=500 , TINY=1 . 0D-20 ) 

INTEGER  i , imax ,  j ,  k 

REAL *8  aamax, dum, sum, vv( NMAX) 

d=l . 

do  12  i=l,n 
aamax=0 . 
do  11  j=l, n 

if  (dabs (a ( i, j )). gt . aamax)  aamax^dabs ( a ( i , j ) ) 

11  continue 

if  (aamax. eq. 0 . )  pause  'singular  matrix  in  ludcmp' 
vv(i)=l ./aamax 

12  continue 

do  19  j=l,n 
do  14  i=l,  j-1 
sum=a( i, j ) 
do  13  k=l, i-1 

sum=sum-a{ i, k) *a(k, j ) 

13  continue 
a( i, j )=sum 
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14  continue 
aamax=0 . 

do  16  i=j,n 
sum=a( i, j ) 
do  15  k=l, j-1 

sum=sum-a( i, k) *a(k, j ) 

15  continue 
a( i, j )=sum 
dum=vv( i) *dabs (sum) 

if  ( dum . ge . aamax )  then 
imax=i 
aamax=dum 
endif 

16  continue 

if  ( j . ne . imax ) then 
do  17  k=l,n 
dum=a(imax, k) 
a(imax/k)=a( j ,k) 
a  ( j  ,  k ) =dum 

17  continue 
d=-d 

vv ( imax ) =vv ( j ) 
endif 

indx( j )=imax 

if ( a ( j  ^  j ) . eq . 0 . ) a ( j , j ) =TINY 
if ( j . ne. n) then 
dum=l . /a( j , j ) 
do  18  i=j+l, n 

a(i,j)=a(i,j) *dum 

18  continue 
endif 

19  continue 
return 
END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  sl#+] . . . 
SUBROUTINE  jacobn(x,yy Tr dfdx, dfdy , n, nmax) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

INTEGER  n , nmax  7 i 

REAL* 8  xy y ( *) , df dx( * ) , dfdy (nmax, nmax) , zl , z2 , z3  7 ely e2 7 e3 , T 
PARAMETER (el=2 6 500 . ,e2=22200 . ,e3=17200. , zl=l . 412997D21 , 
*z2=l . 93216D16 , z3=l . 598361D12 ) 
do  11  i=l,  4 
dfdx(i)=0 . 

11  continue 

dfdy(lr l)=-zl*dexp( -el/T) 
dfdy (1, 2 )=0 . 
dfdy(l7  3 ) —0 . 
dfdy(l, 4 )=0 . 

dfdy (2, l)=zl*dexp( -el/T) 
dfdy (2, 2)=-z2*dexp( -e2/T) 
dfdy (2, 3 )=0 . 
dfdy(2, 4 )=0 . 
dfdy(3, 1)=0 . 

dfdy (3 , 2)=z2*dexp( -e2/T) 

dfdy (3 , 3 )=-2*z3*dexp( -e3/T) *y (3) 

dfdy ( 3 , 4 ) =0 . 

dfdy (4 , 1)=0 . 

dfdy (4 , 2)=0 . 

dfdy (4, 3)=2*z3*dexp(-e3/T)*y(3) 

dfdy(4 , 4 ) =0 . 

return 

END 

SUBROUTINE  derivs (x, y , T, dydx ) 

IMPLICIT  REAL* 8  (A-HyO-Z) 
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REAL*  8  x,y(*)  ,dydx(*) , zl, z2, z3, el, e2,e3,T 
PARAMETER (el=2 6 500 . / 62=22200 . , e3=17200 . , zl=l . 412997D21/ 
*z2=l . 93216D16 , z3=l . 598361D12) 
dydx(l)=-zl*dexp( -el/T)*y(l) 

dydx(2)=zl*dexp( -el/T) *y(l) -z2*dexp( -e2/T) *y ( 2 ) 

dydx(3 )=z2*dexp( -e2/T) *y(2) -z3*dexp( -e3/T)*(y(3 ) **2) 

dydx ( 4 ) =  z  3  *dexp (-e3/T)*(y(3)**2) 

return 

END 

C  (C)  Copr ♦  1986-92  Numerical  Recipes  Software  sl#+] . . . 
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Bijlage  B  Visualisatie  van  computerberekeningen 


Figuur  B.l:  Het  computer  grid  van  de  TNO-PML  Cook-off  buis. 


TIME  COMPLETED  IN  THIS  STEP  i  .  77SB-fOJ  TOTAL  ACCUMULATED  TIMS  1.776B-*0J 
ABACUS  VERSION:  5.4-9  DATE:  19-FBB-9S  TIME:  11:47:55 
STEP  2  INCREMENT  85 


Figuur  B.2:  Berekende  temperatuurcontourplot  van  van  TNO-PML  Cook-off  buis  gevuld 
met  HMX-PBX  vlak  voor  de  thermische  wegloopreactie. 
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Figuur  B.3:  Computermodel  van  een  raketmotor.  De  grijze  kleur  stelt  de  stalen  casing  en 
de  SAD  voor.  Het  gele  de  raket  en  ontstekerstuwstof  Oranje  is  een  HTPB- 
warmteprotectie,  blauw  lucht  groen  een  inhabitor  en  wit  zijn  twee  o-ringen. 
Het  rode  kleine  stukje  vormt  de  ‘ bottle  neck’  voor  het  warmtetransport  naar 
de  ontsteker. 


TIKE  COMPLETED  IN  THIS  STEP  5.6008^03  TOTAL  ACCUMULATED  TIKE  5.600E^03 
ABACUS  VERSION:  5.5-1  DATE:  23-EB3-96  TIME!  15:01:13 
STHP  1  INCREMENT  95 


Figuur  B.4:  Temperatuurcontourplot  van  een  raketmotor  na  5600  sec.  bij  een  opwarm- 
snelheid  van  0,05  °  Cl  sec.  en  Xb=  10  Jim  K  s*. 
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DIKE  I»  THIS  STEP  E.600B-W  TOTU.  ACCUJKJLATBD  TOT  ,.««**> 

ABACUS  VERSION:  S.5-1  DATE i  23-FSB-96  TIKE:  1S,2S,>. 

STEP  1  INCRBMBT1T  9  5 


Figuur  B.5:  Temper atuurcontourplot  van  een  raketmotor  na  5600  sec.  bij  een  opwarm- 
snelheid  van  0,05  0 Cl sec.  en  45  Jim  K  s2. 


TIME  COMPLETED  IN  THIS  STEP  5.S00B-fO3  TOTAL  ACCUMULATED  TIME  5.600B-I-0J 
ABACUS  VBRSICN:  5.5-1  DATE:  23-FEB-96  TIME:  15:20s3< 

STEP  1  INCRBMEflT  9  5 


Figuur  B.6:  Temper  atuurcontourplot  van  een  raketmotor  na  5600  sec.  bij  een  opwarm- 
snelheid  van  0,05  °Clsec  en  \fr=  390  Jim  K  s2. 
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TIMB  COMPLETED  IN  THIS  STEP  1 , 964B-05  TOTAL  ACCUMULATED  TIME  1.9S4B-05 

^ - *  ABACUS  VERSION s  5.5-1  DATEi  23-FBB-96  TIMB;  15:08:18 

STEP  1  1NCSBM&HT  SB 

Figuur  B.7:  Temper atuurcontourplot  van  een  raketmotor  na  196  400  sec.  bij  een  op - 
warmsnelheid  van  3,3  °C/uur  en  Xjy=  10  Jim  K  s^. 


TIME  COMPLETED  IN  THIS  STEP  1.964B-f05  TOTAL  ACCUMULATED  TIME  1.964B-05 
ABACUS  VERSION:  5.5-1  DATE:  2J-FB1-96  TIMB:  10:1S:CB 
STEP  1  INCRBMHtrr  68 


Figuur  B.8:  Temper  atuurcontourplot  van  een  raketmotor  na  196  400  sec.  bij  een  op- 
warmsnelheid  van  3,3  °Cluur  en  X^=  45  Jim  K  s^. 
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TIMB  CCMP LBTED  IN  THIS  STEP  1.964B-f05  TOTAL  ACCUMULATED  TIME  1.9S4B+05 
ABACUS  VERSIONS  5.5-1  DATBi  23-FBB-56  TIME:  15:40:31 
STBP  1  INCREMENT  SB 


Figuur  B.9:  Temperatuurcontourplot  van  een  raketmotor  na  196  400  sec.  bij  een  op - 
warmsnelheid  van  3,3  °Cluur  en  X^=  390  Jim  K  s 2. 


Figuur  B.10:  Temperatuurcontourplot  van  een  gedeelte  van  een  raketmotor  met  reactieki - 
netiek  van  een  HTPBIAP-stuwstofuit  de  literatuur.  De  rode  plek  geeft  de 
plaats  aan  waar  de  stuwstofhet  eerste  ontsteekt . 
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Bijlage  C  Resultaten  van  karakteriseringsmetingen  en 
simulaties 


Figuur  C.l:  Resultaat  van  een  karakteriseringsmeting  van  de  TNO-PML  Cook-off  buis. 


Figuur  C.2:  Resultaat  van  een  karakteriseringsberekening  van  de  TNO-PML  Cook-off 
buis . 
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